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Summary 

This  project  was  originally  proposed  and  funded  for  a  two-year  period.  A  third  year  of 
research  has  recently  been  proposed  and,  while  its  subject  matter  is  a  logical  extension  to  that 
of  the  first  two  years,  it  is  also  a  change  in  direction  (section  III.3).  Therefore  we  present  here 
results  from  the  first  two  years  collectively,  rather  than  the  second  year  alone.  Those  topics  that 
are  exclusive  to  either  year  1  or  year  2  are  identified.  In  some  cases,  however,  the  topics 
described  in  this  report  were  initially  investigated  in  year  1  and  continued  in  year  2. 

During  the  first  two  years  of  the  project,  empirical  studies  of  scale  relationships  in  the 
retrieval  of  sea  ice  lead  statistics  have  been  undertaken,  as  have  modeling  investigations  of 
atmospheric  influences  on  the  satellite  signal.  Additionally,  we  have  developed  statistical  models 
that  describe  the  scaling  properties  of  leads.  The  empirical  studies  have  been  based  primarily 
on  comparisons  within  and  between  Landsat  and  AVHRR  imagery,  while  the  atmospheric  models 
have  been  specific  to  the  AVHRR.  Submarine  sonar  data  have  been  used  in  the  statistical  model 
development.  Specific  accomplishments  to  date  include: 

(1)  The  parameterizations  of  clouds,  haze,  and  atmospheric  chemical  constituents  in  the 
LOWTRAN  7  radiative  transfer  model  have  been  reviewed.  Atmospheric  temperature  and 
humidity  profiles  for  the  arctic  have  been  constructed  from  Soviet  ice  island  data  and  are 
being  used  in  the  model.  (Year  1) 

(2)  Resampling  methods  have  been  tested  on  simulated,  AVHRR,  and  Landsat  images,  as 
have  the  effects  of  digital  enhancements.  (Year  1) 

(3)  Procedures  for  the  retrieval  of  lead  statistics  have  been  developed  and  applied  to  Landsat 
imagery  successively  degraded  to  more  coarse  resolutions. 

(4)  The  relationship  between  "apparent"  lead  widths  measured  along  a  transect  (e.g.,  from 
submarine  sonar  or  as  a  sampling  method  for  satellite  imagery)  and  the  "true"  lead  width 
distribution  has  been  formalized  in  a  statistical  sense,  so  that  one  distribution  may  be 
obtained  from  the  other.  Submarine  sonar  data  have  been  analyzed  in  this  context. 
(Year  1) 

(5)  A  statistical  model  has  been  developed  for  the  retrieval  of  lead  area  fraction  from 
measurements  along  a  line;  e.g.,  a  submarine  sonar  transect  or  a  lineal  sampling  method 
for  satellite  images. 

(6)  The  effect  of  atmosphere/surface  conditions  on  the  AVHRR-measured  radiance  in  the 
thermal  channels  has  been  examined  in  terms  of  thermal  contrast.  Surfaces  include  open 
water,  5  cm,  15  cm,  and  2  m  thick  ice.  Atmospheric  conditions  include  clear  sky  with 
haze,  cirrus,  and  low-level  ice  crystal  plumes. 

(7)  Additional  KRMS  data  has  been  acquired  for  a  time  period  in  which  DMSP  OLS  data  is 
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currently  being  used  in  the  analysis  of  lead  patterns. 

(8)  Information  on  Arctic  aerosol  optical  depths  during  LEADEX  have  been  collected  by  two 
of  the  investigators  on  the  NOAA  P-3  with  sun  photometers.  This  information  will  be 
used  in  future  radiative  transfer  applications. 

(9)  Two  workshops  for  the  satellite  remote  sensing  investigators  of  the  Leads  ARI  have  been 
hosted  by  this  group  in  Boulder.  Recommendations  for  instrumentation  during 
LEADEX '92  were  made. 

One  workshop  report  and  two  referreed  papers  have  been  published,  and  one  paper  has 
been  submitted  for  publication.  Three  additional  papers  are  in  preparation  for  submission  to 
referreed  journals.  Two  graduate  students  have  been  supported  part-time  over  the  course  of  the 
project. 
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PROJECT  OBJECTIVES 

The  goal  of  this  project  is  to  understand  how  sensor  characteristics,  atmospheric 
properties,  and  surface  conditions  influence  the  detection  and  interpretation  of  sea  ice  leads  using 
Advanced  Very  High  Resolution  Radiometer  (AVHRR)  and  other  satellite  data.  We  seek  to 
determine  the  sources  and  magnitudes  of  errors  inherent  in  the  measurements,  how  data  from 
different  sensors  can  be  combined,  and  how  lead  statistics  change  with  the  different  spatial 
resolutions  of  existing  and  future  sensors. 

In  partial  fulfillment  of  these  objectives,  we  have  defined  which  atmospheric  and  surface 
parameters  are  most  critical  for  lead  detection.  Based  on  model  simulations,  we  have  been  able 
to  better  evaluate  the  importance  of  Arctic-specific  model  parameters;  e.g.,  temperature  and 
humidity  profiles  and  aerosols,  for  the  purpose  of  accurately  simulating  sensor  responses.  Sensor 
characteristics  such  as  spectral  response,  field-of-view,  spatial  resolution,  scan  geometry  and  data 
processing  methods  coupled  with  scene  variability  (solar  zenith  angle,  atmospheric  opacity, 
surface  temperature,  snow  cover,  ice  thickness,  size  with  respect  to  sensor  resolution)  determine 
feature  signatures.  The  effects  of  these  parameters  had  to  be  examined  before  lead  signatures 
could  be  evaluated  in  terms  of  lead  width  and  orientation,  particularly  for  features  that  occur  near 
spatial  and  radiometric  limits  of  sensor  resolution. 

The  NOAA  AVHRR  satellite  sensor  provides  daily,  Arctic-wide  coverage  of  ice 
conditions  at  moderate  resolution  and  low  cost.  These  image  sets  contain  information  that  is  of 
primary  concern  to  research  and  operational  interests  in  the  Arctic.  Although  a  variety  of  studies 
have  examined  various  aspects  of  remote  sensing  of  sea  ice,  essentially  no  work  had  previously 
been  done  to  relate  lead  signatures  observed  in  AVHRR  data  to  lead  characteristics.  This  lack 
of  substantive  verification  work  left  key  questions  unanswered  and  posed  significant  research 
problems  relevant  to  current  lead  investigations.  Specifically,  the  following  questions  motivated 
our  research  during  the  first  two  years  of  the  project: 

Lead  and  Surface  Characteristics :  How  does  lead  detection  depend  on  ice  thickness  for 
given  sets  of  sensor  properties,  surface  temperatures,  and  atmospheric  conditions?  Since 
the  temperature  contrast  between  open  water  and  ice  provides  a  means  to  map  leads  using 
thermal  imagery,  to  what  degree  does  this  contrast  affect  the  apparent  width  of  a  lead  as 
observed  in  an  image  and  our  ability  to  detect  it?  Does  lead  orientation  affect  lead 
detection  when  a  wide-angle  scanning  instrument  such  as  AVHRR  is  used  instead  of  a 
nadir-viewing  sensor  such  as  Landsat?  How  accurately  must  surface  temperatures  be 
measured  to  yield  accurate  lead  calculations? 

Atmospheric,  Boundary  Layer,  and  Solar  Zenith  Angle  Effects:  How  do  these  factors 
combine  with  surface  conditions  and  path  length  to  the  sensor  to  determine  the  thresholds 
of  lead  detection?  What  feedbacks  to  die  atmosphere  do  leads  create,  and  how  will  these 
affect  detection;  e.g.,  ice  crystal  plumes  from  open  leads  that  extend  up  to  -  and  in  some 
cases  through  -  the  top  of  the  inversion  layer?  What  are  the  characteristics  of  "typical" 
polar  atmospheres  (i.e.,  water  vapor  content,  temperature  profiles,  cloud  microphysical 
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characteristics),  how  are  they  treated  in  radiative  transfer  models,  and  how  do  these 
factors  affect  remote  sensing? 

General  Sensor  Considerations :  In  what  ways  might  sensor  scan-angle,  sensor 
calibration,  data  gridding,  and  image  enhancements  influence  the  ability  to  detect  leads? 
Are  lead  statistics  derived  from  image  centers  where  spatial  resolution  is  greatest 
comparable  to  those  derived  from  image  limbs  where  resolution  is  poorest? 

While  we  do  not  claim  to  be  able  to  answer  all  of  these  questions  completely,  the  results 
of  the  past  two  years’  research  have  given  us  at  least  partial  answers  to  each.  The  data  collected 
during  LEADEX  will  be  used  in  the  validation  of  those  results. 

To  accomplish  our  objectives,  our  approach  includes  both  modeling  and  empirical  studies 
(Figure  1).  The  radiative  transfer  modeling  is  done  for  the  purpose  of  simulating  the  satellite 
sensor  response  under  a  variety  of  atmospheric  and  surface  conditions.  The  empirical  studies 
include:  a  comparison  of  lead  statistics  determined  in  imagery  of  varying  scales;  e.g.,  AVHRR, 
Land  sat,  KRMS  and  OLS  data;  effects  of  different  resampling  methods  and  digital  image 
enhancements  on  lead  detection  in  AVHRR  and  Landsat  data;  use  of  distributions  derived  from 
the  low  resolution  imagery  to  estimate  characteristics  of  the  distributions  obtained  in  the  high 
resolution  images;  and  the  relationship  between  lead  width  and  spacing  statistics  measured  along 
a  transect  to  the  true  distributions. 

Results  from  this  work  will  be  important  to  the  development  and  application  of  lead 
detection  and  mapping  algorithms  proposed  elsewhere  within  the  Leads  ARI.  For  example,  the 
ability  to  more  accurately  access  lead  width  and  spacing  distributions  from  medium  resolution 
imagery  is  crucial  to  the  evaluation  of  large-scale  heat  flux  estimates.  The  modeling  and 
empirical  approaches  to  quantifying  the  relationships  of  scale  discussed  here  are  a  necessary  first 
step  to  operational  lead  analysis  from  satellite  data. 

Satellite  data  used  in  this  study  includes  Landsat  visible,  AVHRR  visible  and  thermal. 
KRMS  passive  microwave,  and  OLS  visible  and  thermal.  Additionally,  lead  statistics  have  been 
derived  from  submarine  sonar  data.  Radiosonde  temperature  and  humidity  profiles  from  arctic 
ice  islands  are  employed  for  radiative  transfer  studies.  Each  data  type  is  described  in  more  detail 
in  the  appropriate  sections. 

This  report  is  divided  into  three  parts.  The  first  details  those  studies  that  relate  directly 
to  retrieving  lead  statistics  from  satellite  imagery:  image  preprocessing,  sensor  field-of-view,  and 
sampling  methods.  The  second  part  describes  radiative  transfer  studies  of  how  the  atmosphere 
affects  the  satellite  signal  and  how  this  in  turn  might  affect  lead  statistics  derived  from  the 
imagery.  The  third  part  of  the  report  summarizes  the  accomplishments  to  date  and  future  plans. 
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EMPIRICAL 
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Atmos.  Surface  Viewing 


Satellite  Radiances 
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Image  Set 


Figure  1.  Overview  of  the  modeling  and  empirical  approaches  to  the  study  of  lead 
mapping  and  relationships  of  scale  in  satellite  data.  The  radiative  transfer  mode!  on  the 
left  is  used  for  sensitivity  studies  and  to  generate  synthetic  imagery.  This  imagery  as  well 
as  Landsat  data  are  degraded  to  coarser  resolutions.  The  simulated  multi-resolution  data 
sets,  along  with  actual  data  from  different  satellite  sensors,  are  used  to  study  the  effects 
of  measurement  scales  on  derived  lead  statistics. 
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The  empirical  portion  of  this  study  involves  the  comparison  of  AVHRR,  Landsat,  and  - 
to  a  limited  extent  -  KRMS  data.  These  sensors  provide  a  broad  range  of  spectral  and  spatial 
resolutions.  In  this  section,  the  potential  effect  of  image  preprocessing  methods  on  lead  statistics 
is  examined. 

1.1  IMAGE  PREPROCESSING  AND  ITS  EFFECT  ON  LEAD  STATISTICS 

Figures  2  and  3  show  the  effect  on  lead  detection  in  AVHRR  imagery  of  different 
resampling  methods  and  different  digital  number  thresholds  for  "lead  /  no  lead"  mapping.  The 
plots  in  Figure  2  of  digital  numbers  along  transects  in  AVHRR  data  processed  using  different 
resampling  schemes  show  some  relatively  small  shifts  in  digital  number  (DN)  within  leads  and 
adjacent  to  leads,  which  would  yield  a  change  in  estimated  lead-covered  area  (such  as  at  transect 
location  292-293  in  Figure  2a).  For  small  leads,  the  potential  exists  to  perhaps  mask  the  lead 
completely  (as  at  transect  location  279  in  Figure  2a).  Figure  3a  shows  results  using  a  threshold 
that  detects  small  and  large  leads.  The  threshold  used  for  Figure  3b  detects  only  the  larger  leads. 
The  effect  of  the  "smoothing"  interpolations  (bilinear  and  cubic  convolution)  is  to  eliminate  about 
40%  of  the  smallest  leads  using  this  threshold  detection  scheme.  Figure  4  summarizes  the  effect 
of  threshold  and  resampling  scheme  on  total  lead-covered  area  (in  this  case,  the  percent  of  an 
AVHRR  image  covered  by  leads).  The  effect  of  using  different  resampling  schemes  is  small 
compared  to  the  effect  of  choosing  different  thresholds.  Thus,  while  Figure  3  suggests  a 
substantial  reduction  in  small  leads  when  interpolations  are  used  vs.  nearest  neighbor,  the  effect 
of  the  loss  of  these  small  leads  on  total  lead  area  is  relatively  small.  However,  it  is  also  worth 
noting  that,  since  the  turbulent  fluxes  from  a  given  lead-covered  area  increase  as  the  proportion 
of  lead-covered  area  in  narrow  leads  (as  opposed  to  wide  leads)  increases,  any  systematic  shift 
by  the  processing  scheme  toward  a  reduction  in  small  leads  may  need  to  be  considered  when 
calculating  regional  estimates  of  surface  fluxes. 

To  further  examine  the  effects  of  different  interpolation  schemes,  "synthetic"  images 
containing  a  simulated  lead  or  lead  complex  with  different  shapes  and  dimensions  (e.g..  the 
patterns  shown  in  Figure  5)  are  used.  Four  lead  types  were  created  (e.g.,  "Type  1",  etc.  in  Figure 
5):  Type  1  represents  narrow  open-water  leads;  Type  2  wider  open-water  leads;  Type  3  leads  with 
new  and  young  ice  growing  from  the  edges  of  the  leads;  and  Type  4  representing  leads  with  new 
ice  building  up  along  one  side  of  the  leads.  These  images  were  then  resampled  using  bilinear 
interpolation  and  nearest-neighbor  methods.  Resampling  was  applied  to  a  45°  rotation  of  the 
original  images.  Table  1  lists  the  change  in  lead-covered  area  resulting  from  the  resampling  and 
determined  using  two  different  thresholds.  A  DN  of  10  (a  typical  reflectance  for  open  water)  was 
chosen  to  represent  open  water  and  a  DN  of  85  was  used  for  ice-covered  pixels.  The  effects  are 
clearly  dependent  on  the  threshold  chosen  to  define  the  cut-off  between  whether  to  consider  a 
pixel  as  part  of  a  lead  or  not.  Both  the  bilinear  and  nearest  neighbor  interpolations  cause  large 
reductions  in  the  area  assumed  to  be  open  water.  When  a  threshold  is  chosen  at  a  higher 
reflectance,  such  as  would  be  appropriate  to  detect  all  pixels  with  some  open  water  or  thin  ice 
in  them,  then  the  percentage  of  pixels  with  some  lead-covered  area  increases  substantially  when 
bilinear  interpolation  is  used. 
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In  addition  to  affecting  the  calculation  of  lead-covered  area,  the  different  resampling 
schemes  affect  the  appearance  of  leads  in  an  image,  and  thus  the  ability  of  an  interpreter  or 
automated  pattern-recognition  scheme  to  detect  the  leads.  In  particular,  nearest  neighbor 
resampling  tends  to  break  up  the  linearity  of  leads.  Figure  6  shows  a  lead  complex  in  Landsat 
imagery  resampled  using  nearest  neighbor  (Figure  6a),  with  an  attempt  at  reconstituting  the 
linearity  using  a  median  filter  (Figure  6b).  As  with  most  enhancements,  the  improvement  is 
subjective,  but  the  lead  patterns  appear  more  well  defined  in  the  median-filtered  image. 

Additional  insight  into  the  effects  of  spatial  resolution  and  sensor  properties  can  be  gained 
by  comparing  colocated  imagery  from  different  sensors.  Figure  7  shows  registered  Landsat  MSS 
imagery  (top)  and  AVHRR  data  (bottom)  for  a  portion  of  the  Beaufort  Sea.  While  the  same 
general  lead  structure  is  apparent  in  both  images,  the  ability  to  detect  the  smaller  leads  is 
considerably  reduced  in  the  AVHRR  image  (maximum  spatial  resolution  of  1.1  km)  vs.  the  80  m 
Landsat  image.  This  effect  of  spatial  resolution  on  lead  detection,  as  well  as  the  effect  of 
different  spectral  information  on  lead  mapping,  is  also  illustrated  in  Figure  8,  which  shows  a 
subsection  of  the  Landsat  and  AVHRR  images,  as  well  as  a  KRMS  strip  superimposed  on  the 
Landsat  data.  This  comparison  of  how  leads  are  represented  in  visible-band  wavelengths  (the 
Landsat),  thermal  (AVHRR),  and  passive  microwave  (KRMS),  points  out  the  problems  of 
intercomparing  lead  statistics  derived  from  different  sensors.  In  this  example,  threshold  detection 
of  lead-covered  area  in  the  three  data  types  yields  1.1%  lead-covered  area  in  the  Landsat,  12  8% 
in  the  AVHRR  (which  includes  apparent  low  cloud  with  substantially  warmer  temperatures  than 
the  ice  surface),  and  5.4%  lead-covered  area  in  the  KRMS  image. 

These  different  representations  are  perhaps  better  represented  by  comparing  transects 
through  the  imagery.  Figure  9  shows  such  a  transect.  The  lead  located  at  transect  location  32 
is  marked  with  an  arrow  on  the  Landsat  image  in  Figure  8.  The  transect  runs  vertically  through 
the  imagery.  In  this  example,  a  contrast  stretch  was  applied  to  the  AVHRR  data  to  enhance  the 
subtle  DN  differences  in  the  image.  From  these  comparisons  and  a  similar  comparison  of 
transects  in  AVHRR  (unenhanced)  and  Landsat  (Figure  7),  it  is  fairly  clear  that  the  AVHRR  and 
Landsat  reveal  similar  lead  patterns,  but  that  the  number  of  leads  detected,  and  the  image  area 
that  is  considered  to  be  partially  lead-covered,  is  quite  sensitive  to  the  DN  threshold  chosen  to 
define  lead  area.  In  Figure  9,  for  example,  the  large  lead  at  location  32  could  be  defined  as 
having  a  lead  width  from  800  m  to  about  1500  m  in  the  Landsat  data  depending  on  the  DN 
threshold  used,  and  from  300  m  to  1800  m  in  the  AVHRR,  again  depending  on  which  DN  is 
selected.  The  lead  information  contained  in  the  KRMS  data  clearly  is  quite  different  from  that 
shown  in  the  other  image  types. 
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Table  1.  Effect  of  resampling  method  (bilinear  interpolation  =  BL,  nearest  neighbor  =  NN)  and 
digital  number  (DN)  threshold  on  estimation  of  percent  lead-covered  area. 


Fractional  Area  Lead  Coverage  (%) 
(by  Resampling  Method) 

Lead  Type 

Threshold 

None 

Bilinear 

NN 

1 

<  84 

0.65 

1.23 

0.82 

2 

£  84 

1.83 

2.61 

1.77 

3 

<1  84 

1.83 

2.59 

1.77 

4 

£  84 

1.83 

2.60 

1.77 

1 

<  10 

0.65 

0.01 

0.00 

2 

<£  10 

1.83 

1.19 

0.00 

3 

£  10 

1.83 

0.00 

0.00 

4 

<  10 

1.83 

0.64 

0.00 

Transect  Location 


Figure  2.  Change  in  AVHRR  digital  number  (DN)  using  different  resampling  schemes.  Plotted 
data  are  two  transects  (graphs  A  and  B)  through  three  colocated  images  processed  using  nearest 
neighbor  (NN)  (solid  line),  cubic  convolution  (CC)  (dashed  line),  and  bilinear  interpolation  (BI) 
(dot-dash  line). 
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Figure  3.  Comparison  of  number  of  leads  encountered  in  AVHRR  imagery  resampled  using  NN, 
CC,  and  BI  with  two  different  DN  thresholds  (graphs  A  and  B)  used  to  represent  leads. 
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Figure  4.  Change  in  percent  lead  area 
resampling  scheme  and  DN  threshold. 
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Figure  5.  Synthetic  images  used  to  assess  effects  of  resampling  on  different  lead  types. 


Figure  6.  Effect  of  NN  resampling  on  lead  appearance  in  Landsat  imagery  (Fig.  A),  and 
reconstruction  of  lead  appearance  through  the  application  of  a  spatial  filter  (Fig.  B). 
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Figure  7.  Colocated  Landsat  and  AVHRR  imagery  for  the  Beaufort  Sea,  showing  the 
representation  of  the  same  lead  complex  in  the  different  image  types. 


Figure  8.  Colocated  Landsat,  AVHRR,  and  KRMS  +  Landsat  imagery  for  the  Beaufort  Sea.  Tire 
black  arrowhead  in  the  Landsat  image  marks  the  location  of  the  lead  and  transect  presented  in 
Figure  9. 
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Figure  9.  Transect  through  the  imagery  in  Figure  8:  Landsat  (solid  line),  AVHRR  (dashed  line), 
KRMS  (dot-dash  line). 
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1.2  THEORETICAL  ANALYSIS  OF  THE  EFFECT  OF  SENSOR  RESOLUTION 
ON  OBSERVED  FRACTIONAL  AREA  COVERAGE 

While  there  have  been  studies  of  the  effect  of  sensor  resolution  on  parameter  retrieval, 
the  approaches  have  been  empirical  and  have  dealt  only  with  cloud  fraction  (and  a  single  cloud 
type  such  as  cumulus;  e.g.,  Shenk  and  Salomonson  [1972],  Wielicki  and  Welch  [1986])  or  land 
cover  classes  (e.g.,  Woodcock  and  Strahler  [1987],  Townshend  and  Justice  [1988]).  In  the  case 
of  cloud  fraction,  real  and  simulated  data  containing  cloud  fields  were  degraded  in  resolution,  and 
the  fractional  coverage  was  observed  as  a  function  of  scale.  In  the  case  of  land  cover  classes, 
the  variance  of  the  image  was  plotted  as  a  function  of  measurement  scale  for  the  purpose  of 
determining  the  optimal  resolution  for  monitoring.  While  all  of  these  studies  are  useful,  no 
concise  statement  of  the  relationship  between  fractional  coverage  and  sensor  resolution  was  given, 
so  that  the  results  are  difficult  to  generalize  to  other  parameters  such  as  leads.  Figure  10 
provides  an  example  of  visual  changes  in  leads  with  increasing  pixel  size  (and  is  described  in 
more  detail  in  section  1.3).  Here  we  present  an  analytical  approach  to  this  problem  currently 
being  investigated. 

1.2.1  Analytical  Approach 

Our  goal  is  to  determine  the  probability  distribution  for  the  proportion  of  pixels  in  an 
image  that  have  the  characteristic  of  interest;  e.g.,  the  fraction  that  are  cloudy  or  are  sea  ice  leads, 
etc.  This  depends  on  the  distribution  of  the  subpixel  area  fraction  of  the  parameter,  specified  by 
its  shape,  mean  and  variance.  The  variance  depends  on  the  pixel  size  and  the  spatial  structure 
of  the  parameter,  represented  by  its  autocovariance  function.  The  formalization  that  follows  can 
be  applied  to  virtually  any  geophysical  parameter  whose  spatial  distribution  can  be  described  in 
this  manner. 


1.2. 1.1  Indicator  Function 

Let  q(x)  be  a  measurable  property  (e.g.,  temperature)  at  the  point  whose  position  is 
represented  by  the  two-element  location  vector  x,  and  define  £  to  be  any  condition  on  q.  For 
example,  C,  might  be  the  condition  q(x)  <  t- 5,  where  f,  is  the  surface  temperature  and  5  is  some 
threshold  value.  The  indicator  function  l(x)  in  region  R  is  equal  to  1  if  q(x)  satisfies  C,  and  0 
otherwise.  The  fractional  coverage  for  which  q  satifies  £  is  given  by 

P  =  V  l(x)  dx  (1) 


where  AR  is  the  area  of  R  and  is  a  normalizing  factor.  For  the  rest  of  the  discussion  R  is  a 
satellite  image.  The  probability  density  of  I  is/^1)  =  P  and/^0)  =  1  -  P. 

Now  let  q^y)  be  a  measurable  property  of  a  pixel  Z  within  R  centered  at  location  y  (again, 
a  vector).  As  measured  by  the  sensor,  qz  would  be  an  average  over  a  pixel: 
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Qz(y)  =  Az  f  q(x)  dx 

Jziy) 


where  Az  is  the  area  of  the  pixel.  In  reality  the  sensor  point  spread  (transfer)  function  would  play 
a  role  here,  but  this  is  beyond  the  scope  of  this  study.  The  fractional  area  of  R  for  which  q7 
satisfies  £  is  an  estimate  of  the  "true"  fractional  coverage  and  is 

r  -  n-'  £  /*(y)  (2) 

R 


where  N  is  the  number  of  pixels  in  R  and  Iz  is  the  indicator  function  for  the  pixel  based  on  q7, 
defined  in  the  same  way  as  is  I  for  a  point  based  on  q.  Our  goal  is  to  relate  P’  to  P  over  a  range 
of  Az. 


To  determine  P’  analytically  the  probability  density  of  Iz  must  be  known.  It  is  not  trivial, 
and  depends  on  P,  £,  pixel  size,  and  the  way  in  which  objects  satisfying  £  are  distributed  in 
space.  Since  Iz  is  a  function  of  qz,  which  in  turn  depends  on  the  fractional  coverage  within  a 
pixel,  pz,  then  lz  can  be  expressed  in  terms  of  pz.  For  example,  consider  a  cloud  pattern  where 
the  cloud  top  temperature  is  everywhere  the  same.  Let  £  be  a  thresholding  operation  such  that 


lAQz)  - 


if  qz  <  fs-5 
otherwise 


where  ts  is  the  surface  temperature  and  5  is  some  threshold  amount.  (It  is  assumed  that  ts  >  tc, 
the  cloud  top  temperature.)  This  is  equivalent  to 


UPz)  - 


if  ptZ  p 
otherwise 


(3) 


The  expression  pzt  p  states  that  the  fractional  coverage  within  the  pixel  is  greater  than  some 
quantity  p,  which  has  a  value  such  that 

(1  ~P)Tt  +  pTc  <  Trt 


Of  course,  there  is  a  distribution  of  tc,  although  we  do  not  address  this  issue  here.  It  should  also 
be  noted  that  Iz  need  not  have  only  values  of  0  and  1.  It  could,  for  example,  take  on  values  of 
0,  0.5,  and  1  if  £  is  multi-conditional.  So,  based  on  (3),  the  probability  density  of  Iz  is 

f,,(  1)  =  Prob(Pz  Z  p) 

y 0)  =  1  -  Prob(Pz  *  p)  =  Prob(Pz  <  p) 


where  Pz  represents  the  random  variable  subpixel  fractional  coverage  (with  specific  realization 
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pz),  Prob (Pz  ^  p)  represents  the  probability  that  the  fractional  coverage  within  the  pixel  is  greater 
than  some  quantity  p.  Now,  how  is  Prob(Pz  £  p)  determined? 

1.2. 1.2  Distribution  of  Subpixel  Area  Coverage 

For  a  single  pixel  the  fractional  coverage  of  the  geophysical  parameter  is 

Pz(y)  =  Az  (  W  dx 

JAy) 


which  can  be  considered  as  an  estimate  of  P.  After  Stoyan  et  al.  [1989,  p.  184],  the  expected 
value  and  variance  of  Pz  are 

E(PZ)  =  P 

var{Pz)  =  £[PZ-E(PZ)]2 

-  EjlA,1  jz  /(*)  dx  -  P\[AZ 1  fz  /(/)  dx'  -  P] 

■  Az  AcXIIx-x'H)  dx  d/ 


where  k,  is  the  autocovariance  function  for  the  indicator  function. 


If  a  specific  model  distribution  for  Pz  is  assumed,  with  expected  value  and  variance  as 
defined  above,  then  the  density  of  the  pixel  indicator  function  is  also  known.  Here  we  use  the 
Beta  distribution,  a  two-parameter  density  function  defined  over  the  closed  interval  0  <  y  <  1 
often  used  as  a  model  for  proportions: 


=  i 


yo1(  1  -y)**'1 

0, 


r(cc+P) 

mm’ 


a,J3>0;0<y<1 

elsewhere 


The  two  parameters  can  be  determined  from  maximum  likelihood  estimation  based  on  the  mean 
and  variance  of  the  subpixel  fractional  coverage: 

P  -  lE{P;H1-E{PJ)  -var{P,)} 

var\  rz) 


E(Pz)V 

_  [1  -E(PZ)] 


The  shape  of  the  distribution  is  related  to  the  size  of  the  pixel  relative  to  the  spatial 
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structure  (e.g.,  wavelength)  of  the  geophysical  parameter.  In  the  limiting  case  with  very  large 
pixels  relative  to  the  wavelength  of  cloud  elements,  for  example,  most  pixels  would  have  a 
similar  subpixel  fraction  of  cloud  and  the  variance  would  be  very  small.  The  distribution  would 
therefore  have  a  single  peak.  On  the  other  hand,  if  the  pixel  size  is  very  small  then  the 
likelihood  of  pixels  being  either  completely  overcast  or  completely  clear  increases,  the  variance 
increases,  and  two  peaks  are  expected.  This  is  illustrated  in  Figure  1 1  where  the  beta  distribution 
is  shown  for  a  mean  fraction  of  0.2  and  variances  of  0.1,  0.05,  and  0.01. 

Exponential  covariance  is  a  reasonable  model  for  many  geophysical  parameters  and  is 
used  here: 


k^r)  =  p(1  -p)e~ar  (r.a^O)  (4) 

where  a  describes  the  dependence  of  the  covariance  on  the  vector  distance  r=||v-v||.  Implicit 
in  this  expression  is  that  q(x)  is  isotropic.  The  parameter  a  can  be  determined  from  observed 
autocovariances.  Equation  4  can  be  written  in  linear  form  as 

\n[ktf]  =  ln[p(1  -p)]  -  a r 

and  a  determined  through  the  method  of  least  squares  regression.  Examples  of  observed  and 
model  autocovariances  are  given  in  Figure  12. 


1.2. 1.3  Distribution  of  P'  and  its  Error 

As  shown  in  (2),  the  random  variable  P'  is  the  proportion  of  the  N  pixels  in  R  for  which 
/z  =  1 ,  and  its  sampling  distribution  may  be  approximately  normally  distributed  (this  has  not  yet 
been  tested)  with  mean  and  standard  error,  respectively, 

V-  -  4,0) _ 

yi)  y  o) 

N 

if  N  is  large.  The  distribution  function  of  the  fractional  coverage  at  a  particular  sensor  resolution, 
/>,  can  be  determined  from 

FrW  -  X"  ~7=  dx  W 

C\2k 

However,  it  is  not  likely  that  the  Iz  at  one  location  is  independent  of  Iz  at  nearby  locations,  so 
that  the  degrees  of  freedom  must  be  adjusted  to  account  for  spatial  autocorrelation.  This  has  the 
effect  of  increasing  the  variance  of  P'. 
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If  the  true  fractional  coverage  P  is  a  known,  single  value,  then  one  way  of  expressing  the 
error  is  A  =  \P  -  P"  |,  (0  <>  A  £  P),  with  distribution  function 

F„(a)  =  Prob(A£a )  =  W,  dft 


where  the  density  function  fF  is  determined  numerically  from  Ff. 

We  note  here  that  the  preceeding  formalization  of  the  problem  may  not  be  necessary  if 
the  "object"  of  interest  follows  certain  scaling  laws,  such  as  those  used  in  the  definition  of 
fractals.  It  may  be  that  the  relationship  between  changes  in  some  geometrical  descriptor  of  a 
geophysical  phenomenon  and  the  scale  of  measurement  is  adequately  descril  “d  by  the  fractal 
dimension.  This  has  been  shown  to  be  true  of  clouds  [Lovejoy,  1982]  and  land  cover  classes 
[DeCola,  1989]  for  the  relationship  between  area  and  perimeter.  The  stream  length-drainage  area 
relationship  has  also  been  described  in  terms  of  fractals  [e.g.,  Robert  and  Roy,  1990J.  A 
preliminary  investigation  of  lead  length,  L,  at  different  pixel  sizes,  U,  in  the  degraded  Landsat 
data  shown  in  Figure  10  has  also  indicated  a  fractal  characteristic.  This  relationship  is  expressed 
as 


L  -  u'-° 

where  D  is  the  fractal  dimension  (1  <  i  £2).  The  data  in  Table  2  yield  a  fractal  dimension  of 
1.18,  from  which  we  infer  that  lead  lengths  constitute  a  fractal  set. 


Table  2.  Changes  in  open  water  fraction  with  pixel  size. 


Pixel  Size  (m) 

Fraction 

Area  (k 

80 

wgm 

160 

235.8 

320 

196.7 

640 

mmSSmtM 

124.9 

The  above  formalization  is  currently  being  tested  with  a  Poisson  line  process  as  a 
simulated  satellite  image.  The  density  of  the  process  (i.e.,  mean  spacing)  is  1250  m  and  the 
orientations  of  the  lines  (leads)  are  random.  The  lines  are  assigned  thicknesses  (widths) 
following  the  negative  exponential  density  function: 
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where  w  is  lead  width  and  X  is  the  mean  width.  For  the  simulations  X^200m.  One  realization 
of  the  Poisson  line  process  is  shown  in  Figure  13,  where  the  pixel  size  is  30  m.  Figure  14 
illustrates  how  the  scan  angle  effect  can  be  included  in  the  simulated  leads  images. 


Figure  10.  Landsat  MSS  band  4  scene  of  the  ice  pack  north  of  Alaska  in  March  1988.  Area 
covered  in  the  upper  left  image  is  approximately  (80  km)2.  The  remaining  are  degraded  images 
with  pixel  sizes  of  160,  320,  and  640  m. 
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Figure  12.  Autocovariance  of  the  indicator  function  for  two  transects  across  the  Landsat  image 
in  Figure  10.  Shown  are  the  observed  (solid)  and  modeled  (dot-dash)  autocovariances  on  the  left 
and  the  regression  line  fit  to  the  observations  on  the  right. 


Figure  14.  Another  realization  of  a  Poisson  line  process  with  the  AVHRR  scan  angle  effect 
included.  Both  images  are  from  approximately  600  km  off-nadir  on  the  left  side  of  the  image 
and  1350  km  off-nadir  (the  swath  edge)  on  the  right.  In  (a)  all  pixels  are  1.1  km  on  a  side.  In 
(b)  the  scan  angle  effect  is  included. 
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1.3  EMPIRICAL  ANALYSIS  OF  THE  EFFECT  OF  SENSOR  RESOLUTION 
ON  LEAD  WIDTH  DISTRIBUTIONS 

While  there  are  advantages  to  comparing  lead  statistics  derived  from  different  types  of 
imagery  as  shown  in  Figure  8,  such  a  study  would  be  difficult  given  different  acquisition  times, 
wavelength  ranges  of  the  various  sensors,  and  geolocation  problems.  To  alleviate  these  problems, 
comparisons  are  made  across  image  scales  based  on  the  same  initial  image,  where  the  change  in 
resolution  is  obtained  by  modeling  the  transfer  function  between  the  initial  data  and  the  desired 
resolution  and  then  subsampling.  A  spatial  filter  that  estimates  the  point  spread  function  of  the 
Landsat  sensor  is  applied  following  the  methodology  presented  in  Justice  et  al.  [1989].  At  each 
degradation  cycle,  Gaussian  random  noise  is  added  back  into  the  image  to  reduce  the  smoothing 
effects  of  the  filtering  operation.  Images  with  FOVs  of  160,  320,  640,  and  1280  m  were  created 
in  this  manner  (Figure  10).  Data  are  Landsat  Multispectral  Scanner  (MSS)  band  4  (0.5-0.6  pm) 
scenes  of  the  Beaufort  Sea,  March  1988,  with  an  initial  field-of-view  (FOV)  of  80  m. 

Lead  statistics  are  retrieved  from  a  lead/not-lead  binary  image.  To  obtain  the  binary 
image  the  fourth  order  trend  surface  is  removed  from  original  grey  scale  image  [Eppler  and  Full, 
1992]  and  a  thresholding  procedure  is  applied  at  each  degradation  cycle.  One  example  is  shown 
in  Figure  15.  Valid  lead  fragments  are  identified,  where  "valid"  refers  to  a  linear  feature  for 
which  a  meaningful  width  and  orientation  can  be  determined.  Linearity  is  determined  through 
correlation/regression  analysis.  Lead  widths  are  measured  perpendicular  to  the  regression  line, 
at  1  km  intervals,  and  the  slope  of  the  regression  line  is  the  measure  of  the  lead  orientation. 

The  distribution  of  lead  widths  corresponding  to  the  images  in  Figure  10  is  shown  in 
Figure  16.  The  disappearance  of  small  leads  due  to  a  lack  of  contrast  and  the  increase  in  the 
relative  frequency  of  large  leads  as  pixel  size  increases  can  readily  be  seen.  In  this  particular 
image,  a  width  threshold  of  approximately  250  m  divides  those  leads  which  would  disappear  in 
the  320  and  640  m  degradations.  This  change  in  lead  widths  as  a  function  of  pixel  size  affects 
the  total  lead  area  was  previously  illustrated  (Table  2)  for  the  images  in  Figure  10,  where  lead 
area  fraction  decreases  with  increasing  pixel  size.  Orientations  of  leads  can  also  be  expected  to 
change,  if  anisotropy  (i.e.,  a  preferred  orientation)  exists.  An  illustration  of  this  is  shown  in 
Figure  17  for  the  Landsat  image  in  Figure  10.  Results  from  other  Landsat  scenes  show  similar 
patterns  as  do  those  for  the  simulated  images. 

Given  the  changes  in  the  entire  width  distribution  and  in  the  widths  of  individual  leads 
with  increasing  pixel  size  as  illustrated  in  Figure  18,  is  it  possible  to  estimate  the  true  widths  of 
the  leads?  The  criteria  for  how  a  given  lead  will  "grow '  or  disappear  during  image  degradation 
must  include  a  measure  of  its  thermal  contrast  to  the  surrounding  ice.  For  example,  a  narrow, 
open  lead  might  grow  during  the  first  degradation,  accompanied  by  a  drop  in  contrast  due  to 
smoothing,  and  then  disappear  in  the  subsequent  scene.  A  narrow  reffozen  lead,  in  comparison, 
might  disappear  during  the  first  degradation.  An  empirical  transfer  function  that  describes  this 
change  is  currently  being  investigated. 
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Figure  16.  Lead  width  distributions  for  the  images  in  Figure  10.  Widths  are  group- 
bins. 
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Figure  18.  Change  in  lead  widths  of  individual  leads  from  Figure  10  as  a  function  of  pixel  size 
for  those  that  ’•survived"  the  degradation  to  640  m  (top),  and  those  that  did  not  survive  (middle). 
The  bottom  plot  gives  the  means  of  the  top  (solid)  and  middle  (dashed)  plots. 


Part  I:  Image-related  Studies 

1.4  ESTIMATING  WIDTH  DISTRIBUTION  MOMENTS  FROM  PARTIAL  DATA 


35 


Given  that  very  small  features  will  generally  not  be  resolved,  the  issue  then  becomes  one 
concerning  the  possibility  of  using  the  partial  distribution  of  a  particular  statistic  (e.g.,  lead 
widths)  measured  at  low  resolution  to  estimate  the  complete  or  "true"  distribution  of  that 
parameter.  For  example,  assume  that  lead  widths,  x,  follow  a  negative  exponential  distribution 
with  an  unknown  mean,  X.  From  a  sampling  point  of  view  it  is  useful  to  treat  the  distribution 
of  widths  as  discrete  and  address  the  number  n,  of  leads  in  bin  i  that  have  widths  between  x,  and 
x;+w  where  w  is  the  width  of  the  bin: 


*/  = 


Nwo-x,* 

X 


(5) 


where  N  is  the  unknown  total  number  of  leads  in  the  spatial  area.  The  idea  is  that  n,  is  measured 
for  a  few  bins,  and  that  X  and  N  are  estimated.  To  accomplish  this,  (5)  is  rewritten  in  linear  form 
as 


ln(n,)  -  In(^)  -  lx, 

Letting  a  =  In (Nw/X)  and  b  =  X 1  and  solving  for  a  and  b  by  the  method  of  least  squares  with  the 
observed  data,  the  mean  of  the  distribution  and  the  total  number  of  leads  can  then  be  estimated. 

Experiments  with  this  model  show  it  to  be  very  sensitive  to  the  bin  width  and  the  number 
of  bins  in  which  leads  actually  occurred  in  the  low  resolution  imagery.  This  is  not  unexpected 
considering  that  the  entire  range  of  x  is  being  estimated  in  the  least  squares  model  by 
observations  in  only  one  part  of  its  entire  range  (such  extrapolation  is  not  recommended).  The 
problem  can  be  alleviated  somewhat  by  including  a  few  observations  of  small  lead  widths;  for 
example,  measurements  from  a  few  SAR  images  within  the  same  area. 

1.5  LINEAL  METHODS  OF  ESTIMATING  LEAD  PARAMETERS 
1.5.1  Lead  Widths  Measured  Along  a  Transect 

Perhaps  the  largest  source  of  high-resolution  data  that  is  potentially  useful  for  lead 
statistics  is  that  collected  by  submarine  sonar  over  the  last  three  decades.  For  ice  draft 
information  these  data  are  invaluable.  But  can  they  also  be  used  for  statistics  of  lead  geometries; 
i.e.,  lead  widths  and  spacings?  Lead  width  and  spacing  statistics  have  been  examined  in  two 
sonar  data  sets.  Both  are  in  the  Canada  Basin  (Figure  19);  one  in  August  of  1970  and  one  in 
October  of  1978.  Tables  3  and  4  show  the  results.  How  realistic  are  these  data?  If  leads  are 
conceptualized  as  linear  features  of  some  width,  then  crossing  the  lead  at  any  angle  other  than 
perpendicular  to  its  local  orientation  will  result  in  an  overestimate  of  the  actual  width. 
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Table  3.  Lead  widths  (m)  and  standard  deviations  (below)  in  the  QUEENFISH  data  by  region, 
maximum  draft  (cm),  minimum  width  (m).  Also  given  are  the  maximum  lead  widths 
encountered.  No  statistics  are  given  if  fewer  than  20  leads  were  found  in  a  region. 

MAX  DRAFT/  REGION 


MIN  WIDTH 

A 

B 

C 

D 

E 

F 

G 

H 

30/3 

26.5 

33.1 

21.3 

26.5 

47.7 

_ — 

18.5 

23.8 

2.1 

2.9 

2.4 

2.6 

7.1 

2.6 

3.7 

20 

60.5 

3.6 

66.1 

5.8 

54.8 

5.6 

5.2 

100.4 

10.4 

9.6 

50 

200 

Maximum: 

88.2 

5.3 

109.7 

10.8 

228 

502 

130 

283 

327 

70 

90 

70/3 

32.6 

54.8 

27.5 

32.1 

49.4 

36.4 

34.7 

40.0 

3.3 

3.8 

2.0 

2.8 

7.0 

8.3 

3.8 

3.7 

20 

63.7 

78.1 

55.6 

54.2 

73.7 

55.9 

59.0 

60.4 

6.8 

5.2 

3.9 

4.2 

9.9 

13.9 

6.4 

5.4 

50 

106.1 

129.2 

106.5 

104  .2 

129.5 

_ 

98.1 

105.5 

200 

Maximum: 

14.3 

8.5 

7.9 

9.5 

19.0 

10.3 

9.1 

885 

526 

255 

294 

374 

537 

257 

227 

100  /  3 

28.9 

40.8 

32.9 

30.9 

50.5 

32.4 

36.4 

40.0 

2.4 

2.5 

2.5 

2.2 

6.5 

6.6 

3.7 

3.6 

20 

59.3 

70.4 

70.2 

53.6 

74.1 

— 

57.9 

62.8 

5.4 

4.1 

5.5 

3.8 

9.1 

5.9 

5.6 

50 

102.4 

126.6 

124.5 

98.6 

123.0 

98.6 

109.6 

11.9 

7.8 

11.2 

8.1 

15.9 

10.0 

9.8 

Maximum : 

886 

533 

763 

299 

374 

539 

263 

269 
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Table  3,  cont.  Lead  widths  (m)  and  standard  deviations  (below)  in  the 
northeastern  Canada  Basin  data  by  region,  maximum  draft  (cm),  minimum  width 
(m).  Also  given  are  the  maximum  lead  widths  encountered.  No  statistics  are 
given  if  fewer  than  20  leads  were  found  in  a  region. 


MAX  DRAFT/  REGION 


MIN  WIDTH 

I 

J 

K 

L 

M 

30/3 

100.8 

4.9 

2.4 

— 

1.7 

20 

50 

200 

Maximum: 

2181 

153 

83 

42 

70/3 

30.8 

14.4 

17.6 

8.7 

15.9 

16.1 

2.0 

3.3 

1.3 

4.2 

20  132.8  62.5  65.8 

81.4  9.3  13.7 


50 

200 

Maximum: 

2274 

234 

342 

72 

510 

100  /  3 

28.7 

17.9 

18.8 

9.4 

16.5 

12.1 

2.1 

2.7 

1.1 

4.5 

20 

104.6 

61.6 

68.8 

— 

84.7 

50 

200 

52.9 

7.4 

10.8 

29.9 
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Table  4.  Lead  spacings  in  the  QUEENFISH  data,  by  region,  maximum  draft  (cm),  minimum  width  (m). 
S=mean  spacing  (m)  with  standard  deviation  below,  N=number  of  leads  per  kilometer.  Also  given  are 
the  maximum  spacings  encountered.  No  statistics  are  given  if  fewer  than  20  leads  were  found  in  a  region. 


MAX  DRAFT/  REGION 


MIN 

WIDTH 

A 

B 

C 

D 

E 

F 

G 

H 

30 

/  3 

S: 

462.9 

428.5 

929.9 

739.4 

685.9 

— 

2190.9 

2551.0 

39.9 

41.7 

118.9 

82.2 

155.1 

402.5 

759.2 

N: 

2.01 

2.16 

1.06 

1.31 

1.38 

0.45 

0.33 

20 

Ss 

1240.1 

1014.0 

2821.6 

1897.2 

1061.8 

— 

— 

150.1 

125.6 

601.9 

307.5 

240.6 

N! 

0.73 

0.92 

0.31 

0.51 

0.90 

*L  J  «•)  4L  •  v 

318.7 

503.1 

*?■!  ^  J  •  O 

1436.6 

N: 

0.36 

0.41 

0.17 

200 


Maximum : 

4893 

5842 

9321 

8594 

7062 

13523 

24321 

70/3 

S: 

339.7 

369.8 

364.1 

513.0 

504.7 

553.1 

859.6 

737.7 

24.9 

30.2 

34.6 

41.2 

102.8 

108.2 

100.1 

61.9 

N: 

2.67 

2.35 

2.56 

1.83 

1.82 

1.70 

1.10 

1.14 

20 

S: 

783.2 

568.2 

926.6 

1050.0 

816.9 

983.4 

1775.5 

1230.3 

70.6 

50.3 

97.9 

102.5 

173.1 

221.0 

264.0 

131.4 

N: 

1.16 

1.53 

1.02 

0.91 

1.14 

0.98 

0.54 

0.67 

50 

Ss 

1818.8 

1238.0 

2935.7 

3226.0 

1868.1 

3757.9 

3029.1 

190.6 

150.4 

621.9 

577.5 

515.6 

1200.8 

653.3 

N: 

0.50 

0.73 

0.33 

0.29 

0.48 

0.22 

0.27 

200 


Maximum: 

3807 

4  792 

4777 

4488 

6782 

4996 

6221 

3304 

100  /  3 

S: 

238.1 

213.4 

225.5 

416.3 

443.2 

421.6 

837.0 

722.0 

14.8 

13.1 

17.9 

31.5 

88.4 

84.4 

97.6 

72.3 

N: 

3.73 

3.93 

3.86 

2.24 

2.04 

2.20 

1.13 

1.31 

20 

S: 

597.7 

416.5 

581.6 

866.6 

705.2 

1587.7 

1321.9 

46.4 

33.7 

59.9 

78.7 

147.0 

227.5 

213.0 

N: 

1.51 

2.05 

1.54 

1.09 

1.30 

0.60 

0.72 

50 

S: 

1468.9 

1064.2 

1449.1 

2465.7 

1591.7 

- - 

3603.1 

3251.1 

157.0 

126.1 

261.8 

348.7 

397.6 

1158.5 

723.5 

N: 

0.61 

0.83 

0.63 

0.37 

0.60 

0.23 

0.29 

Maximum: 

2897 

2808 

3628 

4478 

6779 

4993 

6218 

8918 
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Table  4,  cont.  Lead  spacings  in  the  northeastern  Canada  Basin  data,  by 
region,  maximum  draft  (cm),  minimum  width  (m).  S=mean  spacing  (m) 
with  standard  deviation  below,  N=number  of  leads  per  kilometer.  Also 
given  are  the  maximum  spacings  encountered.  No  statistics  are  given  if 
fewer  than  20  leads  were  found  in  a  region. 

MAX  DRAFT/  REGION 


MIN 

WIDTH 

I 

J 

K 

L  M 

30 

/  3 

S: 

3381.9 

4567.3 

2477.4 

-  2996.1 

1500.2 

2495.0 

913.9 

1948.1 

N: 

0,15 

0,22 

0.36 

0.17 

20 

— 

— 

— 

- - 

50 

— 

_  _ 

_  _ _ 

200 


Maximum : 

23429 

71342 

42161 

37059 

70/3 

S: 

994.4 

700.1 

928.5 

1763.9 

911.1 

229.6 

129.7 

159.7 

449.3 

169.6 

N: 

0.96 

1.37 

1.01 

0.53 

1.07 

20 

S: 

4476.4 

4587.8 

4479.7 

1396.3 

1513.8 

1513.3 

N: 

0.19 

0.21 

0.20 

50 


200 

— 

— 

— 

— 

— 

Maximum : 

24918 

13508 

10022 

25355 

14562 

100  /  3 

S: 

739.6 

585.6 

635.4 

858.2 

709.5 

130.6 

73.9 

89.7 

197.8 

110.9 

N: 

1.28 

1.65 

1.53 

1.09 

1.36 

20 

S: 

2896.8 

2765.0 

2860.5 

— 

4493.1 

795.8 

623.6 

608.4 

1214.4 

N: 

0.29 

0.35 

0.31 

0.19 

50 


Maximum:  12353  7298  10022  25212  11154 


Figure  19.  Sonar  transect  locations  through  the  central  Canada  Basin  (August  1970)  and 
northeastern  Canada  Basin  (October  1978). 
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A  methodology  has  been  developed  for  estimating  width  distributions  of  linear  features 
from  measurements  along  a  transect  through  a  network  of  such  features.  Both  isotropic  and 
anisotropic  orientations  of  the  linear  features  have  been  considered.  In  the  anisotropic  case,  the 
distribution  of  orientations  of  the  lines  must  be  known.  In  both  cases  if  the  distribution  of  true 
widths  and  the  orientations  are  viewed  as  independent  random  variables  then  the  true  width 
distribution  can  be  determined  from  the  apparent  widths  and  vice  versa.  Furthermore,  if  the  true 
widths  and  the  orientations  vary  jointly  then  the  apparent  width  can  be  determined  from  the  joint 
distribution.  The  geometry  used  in  the  development  of  the  statistical  model  is  shown  in 
Figure  20.  Also  shown  is  an  example  of  the  error  distribution  function  (the  difference  between 
the  true  and  apparent  lead  widths)  for  transects  taken  nearly  parallel  to  the  preferred  orientation 
mean  and  another  perpendicular  to  the  mean  orientation  of  leads  in  the  Landsat  image  of 
Figure  10. 

"True"  lead  width  distributions  in  both  submarine  sonar  and  Landsat  imagery  have  been 
determined  and  apparent  widths  measured  along  transects  were  obtained.  In  the  imagery  it  was 
found  that  lead  orientations  were  not  uniform,  and  width  distributions  generally  followed  the 
negative  exponential  model  although  considerable  variation  was  observed.  The  width 
distributions  measured  from  submarine  sonar  illustrated  that  the  general  shape  is  similar  to  those 
derived  from  overhead  imagery,  but  that  the  errors  in  widths  can  be  significant.  Of  course,  this 
potential  error  is  not  an  issue  if  adequate  two-dimensional  data  are  available.  For  example,  side- 
scan  sonar  will  permit  a  more  accurate  retrieval  of  lead  and  keel  statistics,  although  processing 
such  data  introduces  a  new  set  of  problems.  Sonar  data  with  concurrent  overhead  imagery  from 
aircraft  or  satellite  is  also  a  potentially  valuable  source  of  information.  However,  lead  width 
distributions  derived  from  satellite  data  may  not  be  accurate  since  very  small  leads  are  not 
resolved.  This  is  particularly  true  for  medium  resolution  data  such  as  that  from  the  DMSP  or 
the  AVHRR.  One  solution  might  be  to  retrieve  the  orientation  information  from  satellite  data 
where  small  lead  widths  cannot  be  resolved,  and  the  width  information  from  submarine  sonar 
transects.  Aircraft  using  a  laser  profilometer  would  be  analogous  to  the  sonar  transects. 
Unfortunately,  it  is  not  possible  to  determine  the  actual  error  in  lead  widths  derived  from  much 
of  the  archived  sonar  data.  However,  if  lead  "climatologies"  can  be  compiled  for  various 
locations  and  seasons,  or  if  significant  relationships  between  lead  orientation  and  geostrophic 
winds  can  be  developed,  then  at  least  we  can  determine  the  probable  error. 

Other  applications  of  this  procedure  are  possible.  For  example,  the  formalization  would 
be  the  same  for  lead  and  keel  spacing  distributions  and  their  associated  errors.  An  application 
similar  to  that  of  the  sonar  transect  concerns  sampling  strategies.  Lead  widths  can  be  sampled 
along  a  transect  in  satellite  imagery  or  aerial  photos,  and  this  sample  can  then  be  used  to  infer 
the  characteristics  of  the  true  distribution.  A  similar  application  involves  large-scale  estimates 
of  heat  flux.  Heat  flux  through  leads  is  in  part  a  function  of  fetch,  and  fetch  is  a  function  of  the 
true  lead  width  and  the  crossing  angle  of  the  wind.  If  the  wind  direction  is  constant  as  it  travels 
across  the  network  of  leads,  then  the  distribution  of  fetches  can  be  determined  from  the 
distribution  of  true  lead  widths  as  presented  in  the  previous  sections. 

The  entire  procedure  is  detailed  in  the  Attachment. 


Transect 


Lead 


Error  Distributions 

Transect  Orientations:  2  36  (solid).  0.52  (dashed) 


Error,  m 


Figure  20.  TOP:  Lead  and  transect  geometry  used  to  estimate  the  distribution  of  errors  in  lead 
width  measurements  when  made  along  a  transect;  i.e.,  with  suomarine  sonar.  BOTTOM:  The 
error  distributions  are  shown  for  two  different  transects  across  a  lead  network  that  exhibited  a 
preferred  orientation:  one  oriented  almost  parallel  to  the  mean  lead  orientation  and  the  other 
almost  perpendicular. 
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1.5.2  Lead  Fractional  Area  Coverage  Measured  Along  a  Transect 

In  the  previous  section  sea  ice  leads  were  modeled  as  a  Poisson  line  process  and  a  method 
was  was  presented  for  estimating  the  "actual"  lead  width  distribution  from  a  width  distribution 
measured  along  a  transect  (the  "apparent"  lead  width  distribution).  The  purpose  of  this  section 
is  to  extend  that  work  to  the  estimation  of  the  fractional  area  coverage  of  leads  from 
measurements  along  a  line.  This  line  could  be  a  submarine  transect  under  the  ice  or  on  a  satellite 
image.  Two  methods  from  stochastic  geometry  are  presented  although  applications  have  not  yet 
been  performed.  While  the  problem  is  illustrated  with  respect  to  sea  ice  leads  it  should  be  noted 
that  the  general  lineal  method  described  here  is  applicable  to  any  geophysical  parameter  with 
known  (or  determinable)  spatial  structure. 


1.5.2. 1  Method  for  a  Poisson  Line  Process 

If  in  fact  leads  can  be  accurately  modeled  as  a  Poisson  line  process  then  one  approach 
to  estimating  the  area  fraction  is  to  use  the  lead  width  distribution  [Miles,  1964J: 

ff  =  Prob[0  is  covered\  -  1  -e*w  0) 

where  O  is  an  arbitrary  origin  and  x  is  the  density  of  the  process.  The  overall  mean  line  (lead) 
thickness  (width),  w',  is  defined  as 


W  =  k'  Jjv/(9)d9 


where  w'(0)  is  the  mean  thickness  of  lines  with  orientation  0  (0  £  0  <  n).  This  applies  to  lines 
oriented  isotropically;  i.e.,  with  a  uniform  distribution:  f0(B)=Tt  l  where  fB  is  the  probability 
density  function  for  lead  (line)  orientations.  For  anisotropic  thick  lines  then 

W  =  jV(0)</Fe(8) 


where  dFe(0)  =  /e(0)  d9,  and  Fe  is  the  cummulative  distribution  function  for  orientations.  The 
lead  width  distribution  can  be  determined  as  in  Key  and  Peckham  [1991]. 

1. 5. 2.2  General  Lineal  Method 

The  general  expression  for  the  estimate  of  the  fractional  area  coverage,  /?',  of  a 
geophysical  parameter  whose  actual  fractional  coverage  is  p,  regardless  of  the  the  spatial  structure 
of  that  parameter  is 

P'  «  Hu1  f  /(*)  dx 

V  U 
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The  factor  Pu  is  needed  for  normalization  and  depends  on  the  shape  of  the  structuring  element 
U;  for  example  it  may  be  the  length  of  a  line  or  the  area  of  a  square.  The  indicator  function  I(x) 
for  the  underlying  function  q(x )  at  coordinate  vector  location  x  is  1  if  q(x)  satisfies  some 
condition  and  0  otherwise. 

The  expected  value  of  p'  is  p  and  its  variance  is 


ff  =  E(/f-p)s  =  rffJ^Hx-yDdxdy 


where  k,  is  the  autocovariance  function  of  the  indicator  function  I  [ Stoyan  et  al.,  1989]. 

For  measurements  along  an  array  of  n  parallel  lines  the  corresponding  expressions  are 

ff  -  — 

nl  (2) 


varip)  =  {nf)  z  f  f  kl\\x-y\\)dxdy 

JnlJnl 


where  L  is  the  total  length  of  lines  each  of  length  /. 

Exponential  covariance,  as  described  in  the  previous  section,  is  used  here: 

H1)  =  PC'  -P)®""  (r.a^O)  (3) 

where  a  describes  the  dependence  of  the  covariance  on  the  vector  distance  r=||x-y||.  hnplicit 
in  this  expression  is  that  q(x)  is  isotropic. 

For  a  Poisson  line  process  ct=2t/7t  [Stoyan  et  al.,  1989,  pg.  178;  Miles,  1964]  and  the 
estimation  variance  is  approximately 


vaiipf)  * 


2p(1  ~p)(1 

al 


These  models  are  currently  being  applied  to  real  and  simulated  satellite  data.  Results  will 
be  presented  in  the  next  report. 
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Part  H:  RADIATIVE  TRANSFER  MODELING  STUDIES 
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H.l  RADIATIVE  TRANSFER  STUDIES  IN  LEAD  DETECTABILITY 


II.  1.1  Radiative  Transfer  Modeling  and  Model  Validation 

In  order  to  utilize  radiative  transfer  models  most  effectively,  we  have  reviewed  how  they 
treat  factors  such  as  ice  cloud  morphology,  cloud  optical  thickness,  low  level  inversions, 
boundary  layer  effects,  and  other  aspects  relevant  to  the  polar  regions.  A  particular  concern  is 
that  existing  models,  cloud  properties,  standard  atmospheres,  etc.  have  typically  been  developed 
for  low  and  mid-latitude  applications,  and  may  thus  contain  biases  or  shortcomings  when  applied 
to  polar  regions.  Part  of  this  effort  involves  incorporating  polar  atmospheres  and  cloud  properties 
into  the  models.  Arctic-specific  temperature  and  humidity  profiles  have  been  obtained  and  are 
discussed  in  the  next  section.  Unfortunately,  little  information  on  the  microphysical 
characteristics  of  arctic  clouds  is  available.  Data  that  are  available,  however,  are  being 
incorporated  into  the  models  where  appropriate.  This  is  limited  primarily  to  arctic  stratus 
experiments  during  the  early  1980’s,  and  some  measurements  of  aerosols.  Additional  data  is 
expected  to  be  obtained  during  LEADEX  by  instruments  on  the  NOAA  P-3. 

Work  described  in  this  part  of  the  report  relies  heavily  on  simulating  radiances  measured 
by  the  AVHRR  sensor.  To  simulate  radiances  in  the  AVHRR  thermal  channels,  daily 
temperature  and  humidity  profiles  in  each  season  are  used  with  the  LOWTRAN  7  radiative 
transfer  model  [Kneizys  et  al.,  1988].  Radiances  are  modeled  for  sensor  scan  angles  from  0°  to 
60°  in  10°  increments.  The  appropriate  sensor  response  function  is  applied  to  the  calculated 
radiances,  and  radiances  are  then  converted  to  brightness  temperatures.  Atmosperic  chemical 
composition,  background  tropospheric  and  stratospheric  aerosols  for  the  subarctic  winter  and 
summer  models  are  used,  since  no  such  information  is  available  from  the  ice  islands.  The  optical 
properties  of  Arctic  haze  have  not  been  extensively  measured;  model  calculations  [ Blanchet  and 
List,  1983]  show  that  the  volume  extinction  coefficient  of  Arctic  haze  is  generally  of  the  same 
order  of  magnitude  as  that  of  the  tropospheric  aerosols.  Therefore,  the  use  of  tropospheric 
background  aerosols  is  appropriate. 

In  order  to  test  some  sense  of  the  validity  of  the  radiative  transfer  model,  downwelling 
longwave  irradiances  (fluxes)  computed  with  LOWTRAN  were  compared  to  measurements  at 
South  Pole,  Greenland  and  Denver.  Radiosonde  data  from  the  three  locations  were  used  for  the 
calculations.  To  obtain  irradiance,  E,  from  LOWTRAN  (which  outputs  radiances,  L)  the 
calculation  was  done  for  four  angles  (8  streams)  and  employed  the  weighting  function: 

E  =  71(0.3626838 L,  +  0.3137066 Ls  -  0.222381 0L3  *  0.1012285 L4) 

where  the  subscripts  of  L  refer  to  the  sensor  zenith  angles  79.430°,  58.296°,  37.187°,  and  16.201", 
respectively.  Additionally,  a  single-angle  method  (52.5°)  was  tested  and  found  to  give  fairly 
accurate  results.  The  bandwidth  of  3.5  -  50  pm  used  in  the  calculations  corresponds  with  the 
bandwidth  in  which  a  pyrgeometer  measures. 
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The  LOWTRAN  irradiances  differed  from  the  measurements  by  -4.7%  to  +5.6%  . 
Assuming  5%  accuracy  for  the  pyrgeometer  data  these  results  are  acceptable.  Results  for  clear 
sky  are  listed  in  Table  5. 


Table  5.  Comparison  between  modeled  and  observed  longwave  fluxes  under  clear  sky  conditions. 


Station  and  Date 

Measured 
(W  m'2) 

4  Angles 
(W  m‘2) 

Error  (%) 

52.5° 

(W  m'2) 

Error  (%) 

Denver,  9-30-89 

334 

336 

+0.60 

320 

-4.19 

South  Pole,  12-28-86 

128 

124 

-3.13 

117 

-8.60 

South  Pole,  7-31-86 

71 

75 

+5.63 

71 

0.0 

South  Pole,  1 1-6-86 

107 

102 

-4.67 

96 

-10.28 

I 

Calculations  were  also  done  for  cloudy  sky  conditions.  In  comparisons  of  model  results 
and  observations,  the  observed  cloud  fraction,  Ac ,  must  be  considered.  Here  the  model  is  used 
to  estimate  clear  and  overcast  fluxes,  and  these  are  weighted  by  the  observed  cloud  fraction. 

A  remaining  problem  concerns  the  differences  between  the  microphysical  properties  of 
LOWTRAN’s  cloud  models  and  those  of  the  observed  clouds,  which  are  unknown.  Table  6 
gives  the  results  of  the  cloudy  sky  comparisons. 


Table  6.  Comparison  between  modeled  and  observed  longwave  fluxes  under  cloudy  conditions. 


Station  and  Date 
Cloud  Fraction,  Type 

Measured 
(W  nV2) 

4  Angles 
(W  m'2) 

Error  (%) 

52.5° 

(W  m  2) 

Error  (%) 

Greenland,  7-23-90; 

8/8,  stratocumulus 

309 

301 

-2.59 

301 

-2.59 

Greenland,  7-26-90 

3/8,  stratocumulus 

283 

282 

-0.35 

275 

-2.48 

Greenland,  7-1-90 

1/8,  cirrus 

225 

227 

+0.89 

217 

-3.56 

South  Pole,  7-16-86 

8/8,  cirrus 

157 

156 

-0.64 

156 

-0.64 

South  Pole,  11-26-86 

8/8,  cirrus 

126 

129 

+2.38 

122 

-3.17 
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To  examine  the  effect  of  vertical  temperature  structure  on  upwelling  longwave  radiation, 
radiances  in  the  three  channels  were  estimated  using  arctic  mean  and  subarctic  standard  winter 
and  summer  profiles  (described  below)  with  identical  seasonal  surface  temperatures.  The 
maximum  difference  in  radiances  was  0.05  W  m‘2  sr'1  indicating  that  the  vertical  temperature 
distribution  of  the  relatively  dry  arctic  atmosphere  plays  a  very  small  role  in  the  attenuation  of 
upwelling  longwave  radiation. 


ELI. 2  Ice  Island  Data 

Analyses  are  based  on  atmospheric  temperature  and  humidity  profiles  collected  by 
rawinsonde  from  a  Soviet  ice  island  (NP-26),  located  at  approximately  85°  N  170°  W  during 
1983-1987  (Figure  21).  Generally  two  observations  per  day  were  collected  covering  a  vertical 
range  of  0-25  km.  Profiles  that  have  at  least  10  levels  are  retained  in  the  analysis.  Observations 
include  temperature,  dew  point  depression,  wind  speed,  and  wind  direction.  For  the  years 
1986-87  surface-based  cloud  observations  are  also  available.  These  observations  include  low, 
middle,  and  high  cloud  types,  height  of  the  cloud  base,  and  cloud  fraction. 

Only  clear  sky  profiles  are  of  interest  in  one  of  the  studies  below,  and  since  the  satellite 
thermal  radiances  under  cloudy  conditions  will  reflect  cloud  top  temperature  and  a  significant 
amount  of  cloud  cover  will  affect  the  lower  tropospheric  temperature  structure,  clear  sky 
"seasons"  that  differ  in  their  vertical  temperature  and  humidity  structures  are  then  defined.  The 
seasons  are  determined  objectively  with  a  squared  Euclidean  distance  clustering  algorithm;  the 
variables  are  temperature  and  humidity  at  each  level.  To  reduce  the  degree  of  statistical 
dependence  between  levels,  only  one  measurement  per  kilometer  was  used.  The  resulting  seasons 
are  winter;  October  through  March,  summer:  June  through  August,  and  transition:  April.  May. 
and  September.  The  resulting  mean  seasonal  temperature  profiles  for  clear,  cloudy  (greater  than 
75%  cloud  cover),  and  mixed  conditions  are  shown  in  Figure  22. 


Winter 


Temperature  (K) 


Transition 


Summei 


Temperature  (K) 


Figure  22.  Arctic  winter,  transition,  and  summer  temperature  profiles  under  clear,  cloudy,  and 
mixed  conditions  from  a  Soviet  ice  island  located  in  the  Canada  Basin. 
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n.1.3  Contrast  Studies 

Atmospheric  and  view  angle  effects  are  modeled  under  a  variety  of  conditions  to  simulate 
sensor  response  in  order  to  gain  a  better  understanding  of  how  leads  might  appear  in  the  satellite 
data  under  varying  surface  and  atmospheric  conditions.  With  this  knowledge,  it  may  be  possible 
to  estimate  the  smallest  resolvable  lead  under  a  given  set  of  surface/atmospl  c/sensor  resolution 
and  view  angle  conditions.  Emphasis  is  on  the  AVHRR  thermal  channels  (3:  centered  at 
approximately  3.7  pm,  4:  11  pm,  5:  12  pm). 

LOWTRAN  calculations  using  the  data  and  methodology  described  above  were  made  over 
a  range  of  0.55  micron  optical  depths  for  three  prescribed  surface  temperatures  used  to 
characterize  different  lead  types  (open  or  refrozen)  and  a  fourth  to  represent  the  "background" 
scene  consisting  of  2  m  thick  multiyear  ice.  These  temperatures  were  estimated  using  a  surface 
energy  balance  approach  [Maykut,  1982]  to  characterize  the  open  leads  at  271  K,  leads  covered 
by  5  cm  thick  ice  (257  K)  and  15  cm  thick  ice  (248  K).  The  multiyear  ice  has  a  skin 
temperature  of  235  K  in  all  model  calculations.  The  model  was  run  for  satellite  viewing  angles 
of  0°  (nadir),  20°  and  50°  and  accounted  for  differences  in  directional  emissivity  at  each  central 
wavelength  of  the  AVHRR  channels.  Temperature  and  humidity  profiles  were  based  on  the  mean 
January  ice  island  soundings.  All  simulated  radiances  were  converted  to  brightness  temperatures 
for  ease  of  analysis  and  intercomparisons. 

Examples  of  the  results  are  shown  in  Figure  23  for  three  commonly  observed  conditions 
in  the  Arctic,  ice  crystal  precipitation  and  haze  in  the  boundary  layer  and  high-level  cirrus. 
Figures  23a,  23c,  and  23e  are  composites  of  results  at  0°  viewing  angle  showing  the  behavior  of 
channel  brightness  temperatures  (numbered  3,  4  and  5  on  the  plots)  for  simulations  assuming 
each  of  the  surface  temperatures  defined  above.  Figures  23b,  23d,  and  23f  show  results  for  a 
satellite  angle  of  50°.  A  hypothetical  cirrus  cloud  of  varying  optical  depth  was  inserted  in  the 
lowest  1  km  of  the  atmosphere  to  simulate  the  effects  of  "ice  crystal  precipitation"  which  is  a 
common  feature  affecting  the  radiative  balance  of  the  lower  troposphere  in  the  Arctic. 
Regardless  of  the  underlying  surface  type  the  simulated  brightness  temperatures  rapidly  converge 
to  the  248  K  blackbody  radiating  temperature  of  the  "cloud  top"  level  as  optical  depth  increases. 
Note  that  because  all  three  channels  are  centered  in  window  regions  of  the  spectrum,  surface 
temperatures  are  well  simulated  at  zero  optical  depths.  Maximum  brightness  temperature 
differences  exist  between  the  channels  at  around  optical  depth  0.5  to  0.8.  If  satellite  radiances 
can  be  measured  accurately,  these  differences  may  be  used  to  delineate  leads  when  optically  thin 
ice  crystal  precipitation  is  present  in  the  boundary  layer,  but  brightness  temperature  differences 
and  contrast  between  surface  types  is  rapidly  diminished  as  the  optical  depth  increases  beyond 
about  1.0,  making  it  impossible  to  detect  the  underlying  leads.  This  value  is  not  at  all 
unexpected  as  ice  crystal  precipitation  optical  depths  of  up  to  20  have  been  reported  l Curry  et 
at..  1990], 

Figures  23c, d  show  the  results  assuming  that  a  haze  layer  exists  from  the  surface  to  a 
altitude  of  2  km  (the  boundary  layer  aerosol  model  defined  in  LOWTRAN)  over  the  same  range 
of  optical  depths.  Also  shown  in  Figures  23g,h  are  the  brightness  temperatures  for  aerosols 
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(haze)  during  the  day  when  channel  3  exhibits  a  significant  contribution  from  reflected  solar 
radiation.  Multiple  scattered  solar  radiation  adds  significantly  to  the  upwelling  radiances  as 
optical  depths  of  the  haze  layer  increase  while  the  infrared  channels  do  not  show  any  multiple 
scattering  effects.  All  other  input  parameters  to  LOWTRAN  are  the  same  as  for  Figure  23a.  In 
the  haze  simulations  it  is  obvious  that  the  brightness  temperature  differences  between  channel 
pairs  are  small  when  only  thermal  radiances  are  considered  (nighttime)  over  the  entire  range  of 
optical  depths  while  large  differences  exist  between  channel  3,  and  channels  4  or  5  when  solar 
radiation  is  accounted  for  during  the  day.  Also,  the  contrast  between  surface  types  for  each 
channel  remain  significant  over  the  range  of  visible  optical  depths  shown,  particularly  when  only 
thermal  radiances  are  considered.  Gross  differences  between  Figures  23a, b  and  23c,d  are  due  to 
markedly  different  properties  of  the  haze  and  cirrus  models  defined  for  LOWTRAN. 

The  contrast  between  varying  lead  types  and  the  multiyear  ice  pack  can  be  expressed  by  the 
difference  in  simulated  channel  brightness  temperature  for  a  hypothetical  lead  temperature,  say 
open  water,  5  cm  thick  ice  or  15  cm  thick  ice  and  the  surrounding  ice  pack  at  any  optical  depth, 
normalized  by  the  brightness  temperature  of  the  ice  pack  (here  taken  to  be  2  m  thick). 
Figures  24a-f  show  the  normalized  contrast  for  channel  4  assuming  the  conditions  illustrated  in 
Figures  23a-f.  This  relative  measure  of  contrast  may  provide  a  means  to  detect  leads  using 
AVHRR  data  even  at  night,  though  it  is  apparent  that  if  optically  thick  ice  crystal  precipitation 
exists  that  the  contrast  is  essentially  zero  making  it  impossible  to  detect  leads  using  such  a 
method. 

For  cirrus  clouds  brightness  temperatures  again  converged  at  the  blackbody  radiating 
temperature  at  cloud  top  for  large  optical  depths,  but  the  convergence  was  less  dramatic  than  in 
the  case  of  low  level  ice  crystal  precipitation.  Channel  3  minus  channel  4  brightness  temperature 
differences  at  50°  viewing  angle  (not  shown)  were  significant  in  this  case  reaching  nearly  1 2  K 
over  open  leads  at  a  0.55  micron  optical  depth  of  about  2.6  with  the  difference  being  less  than 
5  K  over  the  pack  ice.  Differences  in  these  signatures  when  incorporated  into  image  analysis 
algorithms  in  conjunction  with  ancillary  information  should  be  very  useful  in  the  detection  of 
leads.  Though  the  normalized  contrast  between  surface  types  degrades  with  increasing  cirrus 
optical  depth,  our  simulations  indicate  that  reasonable  signatures  should  be  detectable  for  a  range 
up  to  about  optical  depth  4.0,  that  is,  for  most  cirrus  conditions  expected  to  be  observed  in  the 
Arctic  region. 

Combining  brightness  temperature  differences  with  contrast  will  improve  our  ability  to 
detect  leads  significantly.  Developing  an  algorithm  to  do  so  presents  a  challenge  for  the  future. 
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Figure  24.  Normalized  contrast  between  an  open  water  lead  and  ice  of  two  thicknesses  for  ice 
crystal  precipitation,  boundary  layer  haze,  and  high-level  cirrus  cloud  at  satellite  view  angles  of 
0°  (a,c,e)  and  50°  (b,d,f).  Conditions  are  mean  January  in  the  central  Arctic. 
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The  ability  to  retrieve  surface  parameters  from  satellite  data  in  the  polar  regions  is 
constrained  by  our  limited  knowledge  of  atmospheric  temperature,  humidity,  and  aerosol  profiles, 
the  microphysical  properties  of  polar  clouds,  and  the  spectral  characteristics  of  the  wide  variety 
of  surface  types  found  there.  In  this  section  we  present  results  in  the  retrieval  of  ice  surface 
temperature  (1ST)  from  the  thermal  channels  of  the  Advanced  Very  High  Resolution  Radiometer 
(AVHRR)  sensors  on-board  the  NOAA  series  satellites. 

Sea  and  land  surface  temperature  (SST  and  LST)  retrieval  algorithms  have  been  developed 
using  the  thermal  infrared  window  portion  of  the  spectrum,  with  the  degree  of  success  dependent 
primarily  upon  the  variability  of  the  surface  and  atmospheric  characteristics.  The  general 
approach  to  estimating  surface  temperature  is  to  relate  satellite  observations  to  surface 
temperature  observations  with  a  regression  model.  Lacking  sufficient  observations,  however, 
satellite  radiances  or  brightness  temperatures  can  be  modeled  by  application  of  the  radiative 
transfer  equation.  This  approach  is  commonly  used  for  SST  retrieval. 

To  our  knowledge,  little  effort  has  been  directed  to  the  retrieval  of  the  sea  ice  surface 
temperature  (1ST)  in  the  arctic,  an  area  where  the  first  effects  of  a  changing  climate  are  expected 
to  be  seen.  The  reason  is  not  one  of  methodology,  but  rather  our  limited  knowledge  of 
atmospheric  temperature,  humidity,  and  aerosol  profiles,  the  microphysical  properties  of  polar 
clouds,  and  the  spectral  characteristics  of  the  wide  variety  of  surface  types  found  there.  We  have 
developed  a  means  to  correct  for  the  atmospheric  attenuation  of  satellite-measured  clear  sky 
brightness  temperatures  used  in  the  retrieval  of  ice  surface  temperature  from  the  split -window 
thermal  channels  of  the  AVHRR  sensors  on-board  three  of  the  NOAA  series  satellites.  These 
corrections  are  specified  for  three  different  "seasons"  and  as  a  function  of  satellite  viewing  angle, 
and  are  expected  to  be  applicable  to  the  perennial  ice  pack  in  the  central  Arctic  Basin 
(Figure  21).  We  do  not  develop  a  completely  new  methodology;  instead  we  modify  a  standard 
procedure  for  use  with  arctic-specific  data.  It  is  assumed  that  a  valid  cloud-clearing  algorithm 
exists  and  that  only  clear  sky  radiances  are  being  examined.  The  cloud  clearing  problem  in  polar 
satellite  data  is  not  trivial,  however.  For  a  review  of  polar  cloud  detection  algorithms,  see  Key 
and  Barry  [1989]  and  Sakellariou  et  al.  [1991]. 

For  the  retrieval  of  1ST  a  multi-channel  algorithm  that  uses  empirical  relationships  to 
correct  for  water  vapor  absorption  is  employed: 

Tic*  "  3  +  bT4  +  cTs  +  d[(T4-Ts)SQcQ] 

where  T4  and  T5  are  the  satellite-measured  brightness  temperatures  (K)  in  the  AVHRR  thermal 
channels  and  8  is  the  sensor  scan  angle.  The  coefficients  are  determined  through  a  least  squares 
regression  procedure,  where  surface  temperatures  are  regressed  against  modeled  brightness 
temperatures. 

AVHRR  thermal  channel  radiances  are  simulated  with  LOWTRAN  as  described 
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previously.  Directional  surface  emissivities  for  snow  are  modeled  [Dozier  and  Warren ,  1982]: 
the  single  scattering  albedo  and  asymmetry  factor  in  the  scattering  phase  function  are  calculated 
from  the  Mie  equations  and  the  directional,  wavelength-dependent  emissivities  are  derived  from 
the  delta-Eddington  approximation  to  the  equation  of  radiative  transfer.  The  directional 
emissivities  are  then  integrated  with  the  response  function  for  channel  i : 

e(X,0)4»/(X)  dk 

£•  m  a 


where  e(X,0)  is  the  emissivity  in  direction  0  at  wavelength  X  and  4>  is  the  sensor  response 
function  which  is  0  outside  of  [X^XJ. 

The  use  of  the  rawinsonde  profiles  in  modeling  the  surface  temperature  requires  an 
additional  step  since  the  first  measurement  in  each  profile  is  the  shelter  temperature,  not  the 
surface  temperature.  Therefore,  the  (unknown)  surface  temperature  for  each  profile  is  assigned 
a  series  of  values  representing  the  range  of  possible  surface  temperatures  for  the  observed 
conditions  during  the  month  to  which  the  profile  belongs.  An  energy  balance  model  [Maykut, 
1982]  is  used  to  determine  these  surface  temperatures,  based  on  the  observed  range  of  shelter 
temperatures  and  wind  speeds  (the  mean  ±1  standard  deviation)  in  the  ice  island  data  for  each 
month. 


The  seasonal  dependence  of  the  coefficients  is  illustrated  in  Table  7,  where  coefficients 
from  each  season  were  applied  to  data  from  every  other  season.  Results  are  shown  for  NOAA  9 
and  indicate  errors  between  0.1  K  for  transition  coefficients  with  winter  data  and  0.6  K  when 
summer  coefficients  are  used  with  winter  data.  Similarly,  the  satellite  dependence  of  the 
coefficients  is  shown  in  Table  8  for  summer  conditions.  On  the  average,  errors  ranging  from  0. 1 
to  1.0  K,  depending  on  season,  can  be  expected  when  applying  coefficients  derived  for  one 
satellite  to  data  from  another,  the  smallest  errors  occurring  between  NOAA  7  and  9  coefficients 
and  data.  Using  SST  coefficients  developed  for  the  North  Atlantic  and  the  Greenland  Sea  area 
to  estimate  1ST  would  result  in  an  underestimate  of  up  to  0.7  K,  largest  in  winter  and  at  scan 
angles  of  40°  and  greater.  While  the  sensor  scan  angle  is  included  explicitly  in  the  correction 
equation,  its  effect  in  the  dry  arctic  atmosphere  is  small,  generally  less  than  0.1  K. 

Surface  temperature  measurements  taken  by  a  PRT-5  thermal  radiometer  during  CEAREX 
in  March  1989  were  compared  to  estimated  ISTs  from  NOAA  11  AVHRR  data.  The  mean  1ST 
for  a  sample  of  four  AVHRR  pixels  was  258.9  K  while  the  mean  PRT-5  temperature  (adjusted 
for  an  emissivity  of  0.998)  of  four  consecutive  measurements  one  kilometer  apart  was  259.04  K. 
Given  the  difficulties  in  comparing  the  two  data  sets  these  results  are  encouraging. 

In  summary,  using  the  split  window  channels  and  scan  angle,  the  rms  error  in  the 
estimated  ice  surface  temperature  is  less  than  0.1  K  in  all  seasons.  Inclusion  of  channel  3  (3.7 
pm)  during  the  winter  decreases  the  rms  error  by  less  than  0.003  K.  Overall,  employing  the  1ST 
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coefficients  results  in  increased  accuracy  of  up  to  0.6  K  over  SST  coefficients  developed  for  the 
North  Atlantic  and  the  Greenland  Sea  areas. 

Additional  details  are  provided  in  the  Attachment. 


TABLE  7.  RMS  error  in  applying  coefficients  (NOAA  9)  developed  for  one  season  (left)  to  data 
from  another  (top). 


Coefficients: 

Winter 

Data  from: 
Summer 

Transition 

Winter 

0 

0.403 

0.128 

Summer 

0.587 

0 

0.342 

Transition 

0.117 

0.219 

0 

TABLE  8.  RMS  eiror  applying  coefficients  (summer)  developed  for  one  satellite  (left) 
to  data  from  another  (top). 


Data  from: 

NOAA  7  NOAA  9  NOAA  1 1 

Coefficients: 


NOAA  7  0  0.272  0.655 

NOAA  9  0.296  0  1.017 


NOAA  11 


0.682 


0.961 


0 
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During  the  first  two  years  of  the  project,  empirical  studies  concerning  scale  relationships 
in  the  retrieval  of  sea  ice  lead  statistics  have  been  undertaken,  as  have  modeling  investigations 
of  atmospheric  influences  on  the  satellite  signal.  Additionally,  we  have  begun  the  development 
of  statistical  models  that  describe  the  scaling  properties  of  leads.  The  empirical  studies  have 
been  based  on  Landsat  imagery,  while  the  atmospheric  models  have  been  specific  to  the  AVHRJR. 
Submarine  sonar  data  have  been  used  in  the  statistical  model  development.  Specific 
accomplishments  to  date  include: 

(1)  The  parameterizations  of  clouds,  haze,  and  atmospheric  chemical  constituents  in  the 
LOWTRAN  7  radiative  transfer  model  have  been  reviewed.  Atmospheric  temperature  and 
humidity  profiles  for  the  arctic  have  been  constructed  from  Soviet  ice  island  data  and  are 
being  used  in  the  model. 

(2)  Resampling  methods  have  been  tested  on  simulated,  AVHRR,  and  Landsat  images,  as 
have  the  effects  of  digital  enhancements. 

(3)  Procedures  for  the  retrieval  of  lead  statistics  have  been  developed  and  applied  to  Landsat 
imagery  successively  degraded  to  more  coarse  resolutions. 

(4)  The  relationship  between  "apparent"  lead  widths  measured  along  a  transect  (e  g.,  from 
submarine  sonar  or  as  a  sampling  method  for  satellite  imagery)  and  the  "true"  lead  width 
distribution  has  been  formalized  in  a  statistical  sense,  so  that  one  distribution  may  be 
obtained  from  the  other.  Submarine  sonar  data  from  the  Canada  Basin  during  summer 
and  fall  has  been  analyzed  in  this  context. 

(5)  A  statistical  model  has  been  developed  for  the  retrieval  of  lead  area  fraction  from 
measurements  along  a  line;  e.g.,  a  submarine  sonar  transect  or  a  lineal  sampling  method 
for  satellite  images. 

(6)  The  effect  of  atmosphere/surface  conditions  on  the  AVHRR-measured  radiance  in  the 
thermal  channels  has  been  examined.  Surfaces  include  open  water,  5  cm,  1 5  cm.  and  2 
m  thick  ice.  Atmospheric  conditions  included  clear  sky  with  haze,  cirrus,  and  low-level 
ice  crystal  plumes. 

(7)  Additional  KRMS  data  has  been  acquired  from  D.  Eppler  (NOARL)  for  a  time  period  in 
which  DMSP  OLS  data  is  currently  being  used  in  the  analysis  of  lead  patterns. 

Two  graduate  students  (E.  Ellefsen  and  S.  Peckham)  were  supported  in  part  by  this 
project.  As  a  consequence  of  their  involvement,  one  paper  has  been  published  and  another  will 
soon  be  submitted  for  publication. 

Additionally,  J.  Key  and  R.  Stone  have  just  completed  sun  photometer  measurements  on 
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the  NOAA  P-3  during  LEADEX.  These  data  will  be  used  to  determine  typical  (during 
LEADEX)  aerosol  properties,  in  particular  optical  depths.  This  information  is  critical  in  the 
radiative  transfer  modeling  discussed  earlier. 

Lastly,  two  workshops  for  the  satellite  remote  sensing  investigators  of  the  Leads  ARI  have 
been  hosted  by  this  group  in  Boulder.  Report  on  the  highlights  of  those  meeting  were  sent  to 
C.  Luther  and  T.  Curtin  late  in  1990.  Recommendations  for  aircraft  instrumentation  during 
LEAD  EX ’92  were  made. 

m. 2  Papers  supported  in  whole  or  in  part  by  N00014-90-J-1840 

(See  Attachments) 

Key,  J.  and  S.  Peckham,  1991.  Probable  errors  in  width  distributions  of  sea  ice  leads  measured 
along  a  transect.  J.  Geophys.  Res.,  96(C10),  18417-18423. 

Key,  J.  and  M.  Haefliger,  1991.  Retrieval  of  ice  surface  temperature,  outgoing  longwave 
radiation,  and  cloud  cover  from  AVHRR  data.  WMO  Report  on  Polar  Radiation  Fluxes 
and  Sea  Ice  Modeling,  WCRP-62,  WMO/TD-No.  442,  C22-C28. 

Key,  J.  and  M.  Haefliger.  Arctic  ice  surface  temperature  retrieval  from  AVHRR  thermal 
channels.  J.  Geophys.  Res.,  in  press  (expected  April  1992). 

Key,  J.,  R.  Stone,  J.  Maslanik,  and  E.  Ellefsen,  1992,  The  Detectability  of  Sea  Ice  Leads  in 
Satellite  Data  as  a  Function  of  Atmospheric  Conditions  and  Measurement  Scale.  Annals 
Glaciol.,  submitted  (March  1992). 

Papers  in  Preparation  Acknowledging  N00014-90-J1840  (tentative  titles;  short  versions  are 
given  in  the  previous  sections): 

•  The  Effect  of  Sensor  Resolution  on  the  Observed  Fractional  Area  Coverage  of 
Geophysical  Variables 

•  Estimation  of  the  Area  Fraction  of  Sea  Ice  Leads  from  Measurements  Along  a  Transect 

•  The  Detectability  of  Winter  Sea  Ice  Leads  in  Thermal  Satellite  Data  Under  Varying 
Atmospheric  Conditions 
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During  LEADEX  a  comprehensive  suite  of  measurements  will  be  taken,  many  of  which 
can  be  used  for  the  validation  of  geophysical  retrieval  algorithms  developed  for  use  with  satellite 
data.  The  measurements  will  be  taken  by  a  variety  of  investigators  and  include  atmospheric  and 
surface  physical  and  radiative  properties,  collected  using  ground-based  and  airborne  sensors. 
Additionally,  a  variety  of  satellite  data  will  be  collected  throughout  the  experiment.  During 
year  3  we  will  obtain  from  the  various  investigators  all  data  pertinent  to  the  retrieval  of  lead 
information  from  satellite  data,  primarily  the  AVHRR,  and  use  these  data  in  the  validation  of  the 
retrieval  algorithms  developed  in  the  first  two  years  of  this  project  based  on  radiative  transfer 
modeling  to  simulate  AVHRR  radiances  as  well  as  empirical  models  of  the  dependence  of  lead 
statistics  (geometrical)  on  measurement  scale. 

In  May  1992,  J.  Key,  J.  Maslanik,  and  two  others  at  CIRES  will  participate  in  the 
Canadian  SIMMS  field  program  in  the  Baffin  Strait/Lancaster  Sound  area.  The  purpose  of  the 
program  is  to  collect  data  that  can  be  used  for  the  validation  of  satellite  retrieval  algorithms,  in 
particular  those  related  to  ERS-1  SAR.  Our  objectives  are  to  collect  data  that  can  be  used  with 
visible  and  thermal  satellite  data,  in  particular  the  AVHRR.  The  data  will  be  an  important 
complement  to  the  measurements  taken  during  LEADEX. 

The  objective  of  year  3  is  to  validate  algorithms  and  empirical  models  developed 
previously  for  the  retrieval  of  lead-related  parameters  from  satellite  data.  This  validation  requires 
a  comprehensive  set  of  surface,  aircraft,  and  satellite  data  that  includes  a  complete  (or  as 
complete  as  possible)  characterization  of  the  surface  radiative  and  geometrical  properties  and  the 
atmospheric  physical  and  radiative  properties.  Specifically,  our  objectives  are  to: 

(1)  Compile  a  complete  satellite  validation  data  set  of  surface,  aircraft,  and  satellite 

observations. 

(2)  Use  this  data  set  to  validate  algorithms  previously  developed  (in  the  first  two  years  of  this 

project)  for  the  retrieval  of 

a.  ice  surface  temperature  (1ST), 

b.  thermal  contrast  between  leads  and  the  surrounding  ice, 

c.  lead  width  distributions  as  a  function  of  measurement  scale,  and 

d.  the  "true"  lead  width  distribution  from  low-resolution  measurements  alone. 

Additionally,  we  will  examine  spectral  reflectance  measurements  of  the  surface  and  relate 
them  to  the  AVHRR  shortwave  channels  for  selected  case  studies. 
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RETRIEVAL  OF  ICE  SURFACE  TEMPERATURE,  OUTGOING  LONGWAVE 
RADIATION.  AND  CLOUD  COVER  FROM  AVHRR  DATA 

J.  Key  and  M.  Haefliger 

Cooperative  Institute  for  Research  in  Environmental  Sciences 
University  of  Colorado 
Boulder,  Colorado  80309-0449 

Introduction 

Satellite  data  for  the  estimation  of  radiative  and  turbulent  heat  fluxes  is  becoming  an 
increasingly  important  tool  in  large-scale  studies  of  climate.  The  ability  to  measure 
parameters  needed  for  flux  estimation  in  the  polar  regions  is  constrained  by  our  limited 
knowledge  of  atmospheric  temperature,  humidity,  and  aerosol  profiles,  the  microphysical 
properties  of  polar  clouds,  and  the  spectral  characteristics  of  the  wide  variety  of  surface  types 
found  there.  In  this  paper  we  present  preliminary  results  of  experiments  in  the  retrieval  of  ice 
surface  temperature  (1ST)  and  outgoing  longwave  radiation  (OLR)  from  the  thermal  channels 
of  the  Advanced  Very  High  Resolution  Radiometer  (AVHRR)  sensors  on-board  the  NOAA 
series  satellites.  The  status  of  polar  cloud  retrieval  from  AVHRR  data  is  also  discussed. 
Inadequacies  in  our  ability  to  measure  essential  parameters  are  identified  and  future  directions 
are  suggested. 

Data  and  methods 

Analyses  are  based  on  atmospheric  temperature  and  humidity  profiles  collected  by 
radiosonde  from  the  Soviet  ice  islands  RVAB  (approximately  80°N  160"W)  during  1979-83 
and  EMIO  (approximately  87°N  175°E)  during  1983-1987.  Generally  two  observations  per 
day  were  collected  covering  a  vertical  range  of  0  to  25  km.  Profiles  that  had  more  than  10 
levels  were  retained  in  the  analyses.  Observations  include  temperature,  dew  point  depression, 
wind  speed,  and  wind  direction.  For  the  years  1986-87  cloud  observations  are  also  available 
for  EMIO.  These  observations  include  low,  middle,  and  high  cloud  types,  height  of  the  cloud 
base,  and  cloud  fraction. 

Mean  summer  and  winter  temperature  profiles  for  clear,  cloudy,  and  mixed  conditions 
are  shown  in  Figures  1,  2,  and  3,  respectively.  Cloud  radiative  forcing  is  evident  in  the 
summer  and  winter  clear  and  cloudy  profiles  where  tropospheric  temperatures  are  higher  under 
cloud  cover  in  the  winter  and  lower  under  summer  cloud  cover.  In  addition,  the  summer 
surface  inversion  disappears  under  cloud  cover.  The  standard  subarctic  winter  temperature 
profile  is  shown  in  Figure  4.  Standard  deviation  at  all  levels  was  computed  for  each 
parameter.  The  maximum  winter  standard  deviations  in  temperature  and  dew  point  depression 
were  6.8  K  and  3.3  K,  respectively  between  0  and  10  km.  In  summer  the  corresponding 
values  are  6.0  K  and  5.8  K.  Interpretation  of  dew  point  depression  variability  requires  care 
since  the  mean  water  vapor  content  is  very  low,  generally  less  than  0.6  gm  m  ’. 

For  1ST  and  OLR  retrieval,  the  individual  profiles  are  used  with  the  LOWTRAN-7 
radiative  transfer  model  to  estimate  radiances  in  the  three  AVHRR  thermal  bands  (channels 
3,  4,  and  5  at  3.7pm,  Upm,  and  12pm,  respectively).  To  avoid  contamination  of  solar 
radiation  in  channel  3,  summer  analyses  are  limited  to  nighttime.  To  examine  the  effect  of 
vertical  temperature  structure  on  upwelling  longwave  radiation,  radiances  in  the  three  channels 
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were  estimated  using  the  arctic  mean  and  subarctic  standard  winter  and  summer  profiles  with 
identical  seasonal  surface  temperatures.  The  maximum  difference  in  radiances  was  0.05  W 
m2  sr 1  indicating  that  the  vertical  temperature  distribution  of  the  relatively  dry  arctic 
atmosphere  plays  a  very  small  role  in  the  attenuatiun  of  upwelling  longwave  radiation. 

ICE  SURFACE  TEMPERATURE  RETRIEVAL 

For  the  retrieval  of  1ST  we  employ  a  multi-channel  SST  algorithm  that  uses  empirical 
relationships  to  correct  for  water  vapor  absorption  (Barton  et  al.,  1989): 

T„  -  a0(8)  ♦  £  ak(Q)Tk  *  aw(d)W  (1) 


where  <7,(0)  are  satellite  zenith  angle-dependent  coefficients,  Tk  are  the  satellite-measured 
brightness  temperatures  in  the  thermal  channels,  and  W  is  the  vertical  water  vapor  column. 
Atmospheric  profiles  under  clear  sky  conditions  for  1986-87  are  used  in  the  analysis. 

Coefficients  are  determined  separately  for  summer  and  winter,  and  for  various 
combinations  of  the  three  thermal  channels.  Additionally,  the  scan  angle  dependence  was 
investigated  for  the  winter  months.  Results  are  given  in  Table  1  for  selected  channel  and 
viewing  angle  combinations. 


Table  1.  Regression  coefficients  for  the  estimation  of  1ST  based  on  AVHRR  thermal 
channels. 


1  Coefficients  j 

R2 

Intercept 

Am 

Am 

AhS 

Summer,  T4,  T3;  0=0° 

0.99 

-2.06 

... 

3.01 

-1.99 

Summer,  T3,  T4,  T3;  9=0° 

1.00 

-4.08 

0.39 

2.65 

-2.02 

Winter,  T4,  T3;  0=0° 

1. 00 

-3.08 

— 

3.06 

-2.05 

Winter,  T„  T4,  T3;  0=0° 

1.00 

-2.21 

-0.14 

269 

-1.54 

Winter,  T3,  T4,  T3;  0=20° 

1.00 

-2.01 

-0.16 

2.63 

-1.45 

Winter,  T3,  T4,  T3;  0=40° 

1.00 

-2.11 

-0.15 

2.71 

-1.54 

The  large  coefficients  of  variation  (rt2)  in  Table  1  indicate  that  the  attenuation  of 
upwelling  radiation  in  the  polar  atmosphere  is  minimal,  especially  in  summer.  Either  of 
channels  4  or  5  alone  could  be  used  in  the  estimation  of  1ST.  Not  correcting  for  atmospheric 
effects  would  result  in  an  1ST  underestimate  of  1. 5-2.1  K  in  channel  4  and  1 .9-3.2  K  in 
channel  5,  summer  or  winter.  The  relationship  between  surface  temperature  and  channel  4 
brightness  temperature  is  shown  in  Figure  5. 

REMAINING  PROBLEMS 

The  primary  problem  in  using  the  ice  island  atmospheric  data  set  concerns  the 
difference  between  "surface"  temperature  as  reported,  and  the  skin  temperature.  The  actual 
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height  of  the  first  measurement  is  unknown,  hut  may  be  a  few  meters  above  the  surface.  The 
difference  between  this  value  and  the  skin  temperature  can  be  positive  or  negative.  We  would 
expect  the  difference  to  be  smallest  over  multiyear  ice  in  winter  and  largest  in  the  summer 
when  the  surface  is  melting.  Until  more  complete  data  are  available  this  will  remain  a 
problem.  In  the  discussions  that  follow,  the  "surface"  temperature  refers  to  the  temperature 
at  the  level  of  the  first  measurement.  A  second  problem  concerns  the  unknown  emissivity 
of  the  surface,  as  well  as  its  bidirectional  emittance  characteristics. 

Future  work  will  include  an  investigation  of  methods  for  estimating  the  skin 
temperature  based  on  the  radiosonde  data.  Such  a  method  may  be  simply  an  extrapolation  of 
the  lapse  rate  back  to  the  surface,  or  may  take  an  energy  balance  approach.  A  new  data  set 
acquired  over  Greenland  is  currently  being  analyzed,  and  should  prove  useful  in  this  regard. 
The  CEAREX  '89  and  M1ZEX  '84  data  sets  will  also  be  useful  for  validation.  Multi-angle 
methods  of  1ST  retrieval,  which  are  based  on  the  differential  absorption  due  to  different 
atmospheric  path  lengths  when  the  same  object  is  seen  from  two  different  view  angles,  will 
be  investigated  when  data  from  the  Along  Track  Scanning  Radiometer  (ATSR)  on-board  the 
ERS-l  become  available. 

Status  of  cloud  parameter  retrieval 

Polar  cloud  detection  algorithms  based  on  satellite  data  have  used  spectral,  spatial,  and 
temporal  information,  yet  no  single  method  has  become  generally  accepted  One  visible  and 
one  thermal  channel  do  not  generally  provide  an  unambiguous  solution  to  the  problem, 
especially  for  low  clouds  over  snow  or  ice  surfaces.  Therefore,  snow  and  ice  cover  data  have 
been  used  as  have  passive  microwave  products.  Some  of  the  basic  ideas  of  the  cloud 
detection  portion  of  the  International  Satellite  Goud  Gimatology  Project  (ISCCP)  algorithm 
(Rossow  et  a!.,  1988)  have  been  adapted  for  use  with  polar  AVHRR  and  SMMR  (Scanning 
Multichannel  Microwave  Radiometer)  data  as  shown  in  Figure  6  (Key  and  Barry,  1989).  In 
test  cases  with  synthetic  and  actual  data,  this  algorithm  has  demonstrated  the  importance  of 
additional  spectral  information  (AVHRR  channel  3  and  the  differences  between  3  and  4)  in 
polar  cloud  identification.  Winter  data  requires  a  separate  algorithm,  in  part  due  to  the  lack 
of  visible  data  and  in  part  due  to  the  differences  in  the  vertical  temperature  structure. 

WTiile  this  algorithm  can  provide  rough  estimates  of  cloud  fraction,  a  more  accurate 
determination  of  this  parameter  as  well  as  the  retrieval  of  cloud  top  height  and  cloud  optical 
thickness  requires  knowledge  of  the  temperature  profile.  Such  information  could  be  derived 
from  TOVS.  although  further  validation  over  polar  regions  is  necessary.  Instead,  given 
sufficiently  small  variability  in  the  vertical  temperature  structure  for  a  homogeneous  region, 
mean  monthly  temperature  and  humidity  profiles  may  provide  enough  information  for  cloud 
detection. 

REMAINING  PROBLEMS 

Unfortunately,  the  multispectral  thresholds  used  in  the  cloud  detection  algorithm  were 
chosen  through  experimentation,  and  lack  a  sound  theoretical  basis.  Additionally,  co-locating 
SMMR  (or  SSM/1)  data  with  the  AVHRR  image  is  computationally  intensive. 

Future  work  includes  the  use  of  TOVS  data  for  the  retrieval  of  temperature  profiles  in 
the  polar  regions,  although  this  will  require  a  significant  validation  effort  (which  is  currently 
in  progress  at  the  Laboratoire  de  Meteorologie  Dynamique,  the  University  of  Washington,  and 
the  University  of  Colorado).  With  .he  possible  exception  of  arctic  stratus,  the  microphysical 
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characteristics  of  polar  clouds  and  haze  have  not  been  studied  adequately,  in  spite  of  their 
importance  in  polar  climates.  The  retrieval  of  sub-pixel  cloud  fraction  using  a  multi-channel 
approach  needs  further  exploration. 

Outgoing  longwave  radiation  from  narrow  band  measurements 

Estimates  of  total  longwave  radiation  flux  are  obtained  by  application  of  a  regression 
equation  to  the  thermal  window  radiances  (Ohring  et  al.,  1984).  This  relationship  takes  the 
form 


Tf  -  +  £  aJ'k  +  bkTl 

»-i 


(2) 


where  Tf  is  the  flux  equivalent  brightness  temperature  (K),  and  n  is  the  number  of  thermal 
channels.  It  may  also  be  useful  to  include  the  difference  between  brightness  temperatures  in 
channels  4  and  5,  which  accounts  for  the  effects  of  water  vapor  loading.  This  form  applies 
only  to  nadir  views,  which  may  be  estimated  by  means  of  a  limb-darkening  function. 
Equation  (2)  shows  a  quadratic  form  of  the  relationship,  although  the  linear  form  may  be 
sufficient  and  is  used  here.  An  example  of  this  relationship  is  shown  in  Figure  7,  which 
relates  7}  to  T4. 

As  with  1ST  retrieval,  coefficients  are  determined  separately  for  summer  and  winter. 
Results  are  given  in  Table  2  for  channels  4  and  5,  and  for  all  three  thermal  channels.  Scan 
angle  dependencies  have  not  yet  been  investigated  for  cloudy  conditions. 


Table  2.  Regression  coefficients  for  the  estimation  of  OLR  based  on  AVHRR  thermal 
channels. 


|  Coefficients  j 

R2 

Intercept 

Ah3 

Am 

Aw 

|  Winter,  TJ(  T3 

0.87 

74.33 

1 

t 

1 

1.42 

-0.76 

j  Winter.  T,.  T4,  T, 

1.88 

67.71 

0.32 

066 

-0.29 

REMAINING  PROBLEMS 

While  the  results  presented  above  for  the  three  thermal  channels  are  encouraging,  the 
estimation  procedure  could,  in  theory,  be  improved  by  sampling  the  longwave  portion  of  the 
spectrum  at  additional  wavelengths.  Using  other  thermal  channels  from  sensors  such  as  the 
H1RS  sounder  has  been  investigated  elsewhere  (cf.,  Eilingson  et  al.,  1990)  and  has 
demonstrated  increased  accuracy. 

Future  work  includes  an  investigation  of  scan  angle  dependence  and  limb  darkening 
functions,  and  a  comparison  of  the  modeled  broadband  irradiances  to  those  measured  by  the 
Earth  Radiation  Budget  (ERB)  satellite. 
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Probable  Errors  in  Width  Distributions  of  Sea  Ice  Leads 
Measured  Along  a  Transect 

j.  Key 
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S.  PECKHAM 
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The  degree  of  error  expected  in  the  measurement  of  widths  of  sea  ice  leads  along  a  single 
transect  are  examined  in  a  probabilistic  sense  under  assumed  orientation  and  width  distributions, 
where  both  isotropic  and  anisotropic  lead  orientations  are  examined.  Methods  are  developed  for 
estimating  the  distribution  of  "actual"  widths  (measured  perpendicular  to  the  local  lead 
orientation)  knowing  the  "apparent”  width  distribution  (measured  along  the  transect),  and  vice 
veraa.  The  distribution  of  errors,  defined  as  the  difference  between  the  actual  and  apparent  lead 
width,  can  be  estimated  from  the  two  width  distributions,  and  all  moments  of  this  distribution  can 
be  determined.  The  problem  is  illustrated  with  Landsat  imager;  and  the  procedure  is  applied  to 
a  submarine  sonar  transect.  Results  are  determined  for  a  range  of  geometries,  and  indicate  the 
importance  of  orientation  information  if  data  sampled  along  a  transect  are  to  be  used  for  the 
description  of  lead  geometries.  While  the  application  here  is  to  sea  ice  leads,  the  methodology  can 
be  applied  to  measurements  of  any  linear  feature. 


1.  Introduction 

Polar  sea  ice  is  an  important  factor  in  the  complex 
interaction  of  ocean  and  atmosphere.  Reduction  in  the 
extent  and  thickness  of  sea  ice  due  to  global  warming,  and 
the  consequent  increase  in  the  number  of  cracks  in  the  ice 
(hereafter  "leads"),  is  expected  to  further  increase  global 
temperatures.  This  positive  feedback  is  a  result  of 
reduced  albedo  and  the  increase  in  heat  transfer  from  the 
ocean  to  the  atmosphere.  Model  estimates  indicate  that 
an  increase  of  4%  in  the  area  covered  by  leads  during 
winter  could  produce  a  hemisphere-wide  warming  of  1 
degree  Kelvin  [Ledley,  1988].  Understanding  lead  forma¬ 
tion  processes  as  well  as  the  geographical  and  temporal 
distribution  of  lead  networks  is  therefore  important  to 
studies  of  global  climate. 

Measurements  of  ice  thickness  and  lead  coverage  are 
commonly  made  along  transects  using  upward-looking 
submarine  sonar.  The  footprint  of  the  sonar  beam  used  by 
submarines  to  observe  ice  conditions  is  of  sufficiently  high 
resolution,  often  3  m  or  less  as  opposed  to  25  m  or  more 
for  satellite  sensors,  to  make  it  a  potentially  useful 
instrument  for  gathering  lead  statistics.  In  fact,  sonar 
data  have  shown  that  the  largest  number  of  lead  widths 
are  in  the  10-20  m  range.  Lead  width  statistics  derived 
from  submarine  sonar  data  have  been  reported  in  the 
literature  [e.g McLaren,  1989;  Wadhams,  1981;  Wad  hams 
and  Home,  1980],  but  when  lead  orientation  and  width 
statistics  are  not  available,  the  error  in  these  analyses 
cannot  be  accessed.  This  error  may  be  relatively  small  if 
leads  are  narrow  or  randomly  oriented,  but  may  be  large 
in  the  case  of  wide  leads  or  leads  with  a  preferred  orienta¬ 
tion.  While  it  may  not  be  possible  to  determire  the  actual 
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error  for  much  of  the  archived  sonar  data,  an  important 
question  concerns  whether  or  not  the  data  can  be  used  for 
lead  geometry  statistics,  and  if  so,  what  is  the  magnitude 
of  the  maximum  error.  In  this  paper  a  probabilistic 
determination  of  this  error  is  described,  providing  a 
starting  point  for  the  application  of  stochastic  geometry 
theorems  in  the  analysis  of  lead  geometries.  Errors  in 
statistics  derived  for  other  lead  and  keel  features  are 
discussed  briefly.  While  the  application  is  to  sea  ice  leads 
and  sonar  data,  the  methods  also  apply  to  the  general 
problem  of  sampling  linear  features  along  a  transect. 

2.  Definitions,  Notation,  and  an 
Illustration 

In  the  following  discussions,  notation  follows  that  used 
in  probability  theory,  where  F^z)  denotes  the  distribution 
function  (df)  for  the  population  random  variable  Z  with 
specific  instance  z  (i.e.,  F^lz)  m  P [Z  <  *])  and  /jjte)  is  the 
probability  density  function  (pdf).  Additionally,  E[Z]  and 
VartZ]  are  the  expected  value  and  variance  of  Z. 

The  problem  is  to  relate  a  distribution  of  lead  widths 
taken  along  a  line  perpendicular  to  the  local  orientation  of 
a  lead  (the  ’  actual"  width)  to  the  lead  widths  measured 
along  a  transect  (the  "apparent"  lead  width),  taking  into 
account  lead  orientations,  and  lead  crossing  angles.  As 
illustrated  in  Figure  1,  the  following  continuous  random 
variables  are  defined:  X  is  actual  lead  width,  X  is  appar¬ 
ent  lead  width  measured  along  a  transect,  6  is  lead 
orientation  (0  S  0  <  x),  and  A  is  lead  intersection  angle  (0 
S  A  <  n),  with  specific  realixations  x,  x’ ,  9,  and  a.  Addi¬ 
tionally,  let  4  be  the  transect  orientation  (0  £  4>  <  *)•  The 
position  and  orientation  of  a  lead  within  the  plane  are 
uniquely  specified  by  the  length  of  the  perpendicular  that 
connects  the  lead  to  the  origin,  and  the  angle  that  it 
makes  with  a  fixed  reference  line.  The  intersection  angle 
A  is  measured  between  the  transect  and  the  lead,  anti¬ 
clockwise,  and  is  the  difference  between  their  orientations. 
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Pig.  1.  The  geometrical  relationships  between  a  lead  and  a 
transect.  See  text  for  definition  of  angles  and  length  variables. 


Finally,  define  A'  -  |n/2  -  A  |  to  be  the  crossing  angle  (a' 
in  Figure  1)  measured  between  the  transect  and  a  perpen¬ 
dicular  to  the  lead  orientation  (0  £  A'  £  nJ2).  The  relation¬ 
ship  between  apparent  and  actual  lead  widths  is 


x  = 


(1) 


X 

coeCA') 

where  X  <,  X.  Rearranging  terms,  a  lead  crossing  angle 
can  be  determined  from  the  lead  widths  by 


A’ 


Fig.  2.  Landsat  MSS  band  4  scene  of  the  ica  pack  north  of  Alaska 
in  March  1988.  Area  covered  is  approximately  (80  km)2.  Field-of- 
view  is  80  m 
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Fig.  3.  Lead  width  distribution  for  the  scene  in  Figure  2.  Widths 
are  measured  along  a  perpendicular  to  the  local  orientation  of  the 
lead,  and  are  grouped  in  100  m  bins.  The  mean  width  is  348  m, 
the  standard  deviation  201  m,  and  maximum  width  1376  m. 


The  potential  inaccuracies  of  measuring  lead  widths  along 
a  transect  can  be  illustrated  by  randomly  choosing  a 
transect  orientation  and  location  on  a  satellite  image. 
Here  we  provide  an  example  with  a  Landsat  Multispectral 
Scanner  (MSS)  band  4  (0.5-0.6  pm)  scene  of  the  Beaufort 
Sea,  March  1988  (Figure  2).  The  pixel  size  is  80  m;  image 
size  is  80  x  80  km,  a  subset  of  a  Landsat  scene.  To 
increase  the  sample  size  of  lead  widths  measured  along 
the  transect,  multiple  transects  of  the  same  orientation 
are  placed  randomly  on  the  image.  It  is  assumed  that  the 
pattern  of  leads  is  similar  beyond  the  image  boundaries. 
Processing  of  the  Landsat  data  for  the  retrieval  of  lead 
statistics  is  as  follows.  A  dynamic  threshold  procedure  is 
applied  that  estimates  the  probability  density  function  of 
a  mixture  population  (lead/ice)  for  small  regions  within  the 
image,  and  a  binary  image  results.  Valid  lead  fragments 
are  identified,  where  "valid"  refers  to  a  linear  feature  for 
which  a  meaningful  width  and  orientation  can  be  deter¬ 
mined.  Linearity  is  determined  through  correlation/regres¬ 
sion  analysis.  Lead  widths  are  measured  perpendicular  to 
the  regression  line,  every  kilometer  along  the  lead  length, 
and  the  slope  of  the  regression  line  is  a  measure  of  the 
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Fig.  4.  Lead  orientations  for  the  scene  in  Figure  2.  The  mean 
orientation  is  0.67r  (38.4")  with  standard  deviation  0.87T  (49.8°). 
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Pig.  6.  Lead  widths  from  a  randomly  choeen  transect  across 
Figure  2.  Transect  orientation  is  3.0*  (172°)  or  approximately 
aouth-eo'  .thwest  to  north-northeast  where  the  top  of  the  image  is 
north  rhe  mean  width  is  368  m,  the  standard  deviation  474  m, 
am1  maximum  width  2818  m. 


lead  orientation.  Further  details  of  this  procedure  are 
given  in  Key  et  al.  [1990]. 

The  distribution  of  these  actual  lead  widths  x  is  shown 
in  Figure  3  and  orientations  6  in  Figure  4.  The  mean  lead 
orientation  is  0.67  radians  (38°,  approximately  southeast 
to  northwest  where  the  top  of  the  image  is  north).  For  a 
transect  orientation  $=3.0  radians  (172°,  south-southwest 
to  north-northeast),  the  distribution  of  apparent  lead 
widths  x'  is  illustrated  in  Figure  6,  with  crossing  angles  a' 
shown  in  Figure  6.  The  mean  actual  width  is  348  m  with 
a  standard  deviation  of  201  m,  while  the  mean  apparent 
width  is  368  m  with  a  standard  deviation  of  474  m. 
Additionally,  the  maximum  actual  lead  width  in  the  image 
is  1376  m,  while  the  maximum  width  measured  along  the 
transect  is  2818  m.  With  a  transect  orientation  of  0.13 
radians  (7.4°)  the  difference  between  the  actual  and 
apparent  mean  widths  is  139  m  and  the  maximum  width 
is  2670  m.  From  this  example  it  is  clear  that  significant 
errors  can  result  from  sampling  along  a  transect.  The 
following  section  presents  a  method  to  assess  this  error. 


Fig.  6  Lead  creasing  angles  for  a  transact  a  cross  Figure  2. 
Transact  orientation  is  the  same  aa  in  Figure  5.  The  mean 
creasing  angle  is  1.09*  (62®)  with  standard  deviation  0.23T  (13°). 


3.  Probability  Models 

Theorems  of  stochastic  geometry  that  are  applicable 
here  have  been  developed  through  the  study  of  fibers  as  a 
stationary  random  process  in  the  plane.  If  we  use  this 
analogy  with  lead  networks,  then  after  Stoyan  et  al.  [1987, 
p.  240]  the  df  of  intersection  angles  is 

f**°  sin  (0'  -  $)  dFe(&) 

Fa(  a)  =  * _  (2) 

JJ1  |  sin  (e'-$)|  dFem 

where  FA(a)  is  the  probability  of  intersection  angles 
between  0  and  a,  Fq(9)  is  distribution  function  for  lead 
orientations,  dFe(Q)  =  /q(0)  d0,  a  increases  in  an  anti¬ 
clockwise  sense,  and  in  the  integral  F^fn  +  a)  =  1  +  Fe(a). 
The  pdf  fe(Q)  may  be  an  assumed  mathematical  distribu¬ 
tion  or  may  be  based  on  an  observed  rose  of  direction. 

If  the  leads  are  isotropic  then  the  corresponding 
orientations  have  a  simple  uniform  probability  distribution 
in  the  interval  0  5  0  <  n;  i.e.,  /e(0)=l/n  for  all  0.  In  this 
case  the  distribution  of  intersection  angles  is  independent 
of  the  transect  orientation.  The  probability  of  crossing  a 
lead  that  is  oriented  across  the  transect  (a  — »  n/2)  iB 
greater  than  for  one  running  more  parallel  (a  — »  0  or  n). 
The  associated  intersection  angles  have  density 

fA(a)  ~  W  sin  a  ,  0<a<n 

which  is  not  uniform  but  is  symmetrical  about  a  =  n/2. 
The  corresponding  distribution  function  is 

FA(a)  =  J^a  fA(\y)  dy  =  Vi  (1  -  cos  a)  ,  0Sa<n 

which  is  a  special  case  of  (2)  for  /^(Q)  =  1 In.  In  the 
anisotropic  (preferred  orientation)  case  we  use  (2)  for  the 
distribution  of  intersection  angles,  and  the  corresponding 
densities  are  determined  numerically. 

Two  different  intersection  angles  correspond  to  each 
crossing  angle  so  that  the  distribution  and  density  func¬ 
tions  for  the  crossing  angle  are 

Fa<c O  a  F\A'  £  u'] 

=  PI  |* -A  |  <!  al 

=  P[-a'  <,  (A-l)  S  ''] 

2 

«  Pl(^-a')  SA  <;  i*+a')] 

-  FA** o')  -  PA(|-a'l 

and 

fAda’)  =  =  fA(**al  *  fA(*- a') 

aa  4  4 

which  in  the  isotropic  case  yields  FA\a')  =  sin  a'  and 
fAAa')  =  cos  a’. 

3.1.  Width  Distributions 

An  expression  for  the  joint  density  function  of  the 
apparent  and  actual  lead  widths  can  be  derived.  Suppose 
that  the  joint  pdf  of  X  and  A’  are  known,  which  will  be 
fy. |-(*,a')  =  if  the  two  variables  are  indepen¬ 

dent.  Then  ifV  and  Z  are  two  new  random  variables  that 
are  functions  of  X  and  A'  such  that  Y  =  X  and  Z  =  X  = 
X/cos  A',  then  the  joint  pdf  of  Y  and  Z  can  be  computed 
using  a  standard  theorem  [Poss,  1984,  p.  217]: 
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■  fxx^x*') 

-  /ar,A'U,co8‘1(ij.)]  4  K*')2  -  x2]'^ 

XX 

=  fx<*)  fA'i C08-1(4>]  4  [(X)2  -  xzr*(3) 
x  x 

The  first  of  these  expressions  is  valid  whether  or  not  X 
and  A'  are  independent.  If,  however,  future  research 
indicates  that  large  leads  are  oriented  differently  than 
small  leads,  for  example,  then  the  joint  density  function 
must  be  determined  in  another  manner  (possibly  from 
observations).  Using  (3),  the  pdf  of  apparent  lead  widths 
can  be  obtained: 

&■<*'>  =  fxjC^’X')  dx 

-  J £>*<*>  l/A'tcos-^i,)]  4 

*  K*')2  -x2r»\Xdx  x  <4) 

The  df  of  apparent  lead  widths  can  be  obtained  by  inte¬ 
grating  (4)  or  by  conditioning  on  the  value  of  X,  again 
assuming  that  X  and  A'  are  independent.  The  latter 
method  yields 


numbers  that  must  be  evaluated  at  log  x,  which  becomes 
a  sampling  problem  with  discrete  data;  i.e.,  a  very  large 
number  of  observations  sampled  at  small  intervals  would 
be  necessary  for  any  reasonable  degree  of  success  with  this 
approach. 

Fortunately,  there  is  another  way  to  approach  the 
problem.  Referring  back  to  (4)  and  (5),  fx  and  Fx  can  be 
viewed  as  the  result  of  applying  an  integral  operator  to  fa. 
The  functions  can  be  discretized  as  arrays  and  the  integral 
in  (5)  approximated  as  a  sum: 

J 

Fxii)  =  A  5Z  FA>(j,i)  fx(i),  (9) 

»=1 

i,je[ljf],  x~iA,  y-jA 

where  FAJj,i)  =  F^cos'^i /j)],  N  is  the  number  of  discrete 
observations  and  A  is  the  increment  between  observations. 
If  these  functions  are  expressed  as  matrices,  (9)  becomes 

Fy  =  A  Fa  fx 

fx*±  F  a  Fr 

whose  derivation  is  given  in  the  appendix. 


Fy(x')  =  jf'  Fa-[cob-1(4)]  fjfic)  dx  (6) 

which  is  based  on  the  df  rather  than  the  pdf  of  A'. 

Determining  the  distribution  of  actual  lead  widths  given 
the  apparent  lead  width  distribution  must  be  approached 
differently.  Letting  Y  =  l/cos  A',  then  (1)  can  be  rewritten 

as 

X ’  -  XY 

log  x-  =  log  x  *  log  y 

X  and  A'  are  assumed  to  be  independent,  hence  so  are 
log*  and  log  A',  so  that  log  X'  is  a  sum  of  two 
independent  random  variables.  This  allows  us  to  write 

hog  X  3  hog  X  •  hog  Y  (6) 

or 

*log  X  =  *log  X  *  hog  Y  (7) 

where  the  asterisk  represents  convolution  [Foss,  1984, 
p.  202].  The  Laplace  transform  may  be  used  with  either 
(6)  or  (7).  The  Fourier  transform  may  be  used  with  (6)  but 
not  with  (7)  because  the  pdfs  are  absolutely  integrable  but 
the  df  is  not.  Using  <p (/]  to  denote  the  Fourier  transform  of 
f,  then  with  (6) 


<rt/iog  =  <pt/iog  x’hog  rl  =  ^t/iog  xWog  y] 


hogX  -  9 


9l/iog  xl 


from  which  the  actual  lead  width  distribution  is 


fxfx)  =  1  9'1 

X 


Wiog  xl 
9(/iog  yJ 


(log  x) 


(8) 


The  derivation  of  (8)  is  given  in  the  appendix.  This 
expression  shows  that  when  the  appropriate  transforms 
exist,  then  the  pdf  ofX  is  uniquely  determined  in  terms  of 
the  pdfs  of  X"  and  A'.  However,  it  is  only  useful  for 
computations  if  the  Fourier  transforms  and  inverse 
transforms  shown  can  be  performed  analytically.  This  is 
because  in  practice  the  inverse  transform  is  an  array  of 


3.2.  Error  Distribution 

In  thiB  study  the  error  in  measured  lead  width  is 
defined  as  X-X  (which  is  always  positive),  although  other 
definitions  such  as.XXX'  would  also  be  useful.  Equation  (3) 
allows  us  to  compute  the  distribution  of  the  error  as 
follows: 


Fy.^a)  =  P[X-X  £  a] 

=  X**  fx  *  fojT<x’x')  d*’  <**  >  a  *  0 


For  the  isotropic  case 

Fy.jjfa)  =  fxix)  JZ 
All  moments  of  X-X  can  be  computed  from  /y.y 


/aT 


dx 


4.  Application 

These  models  are  now  applied.  First,  a  lead  width 
distribution  measured  from  sonar  data  is  used  to  estimate 
the  actual  lead  width  distribution,  for  both  isotropic  and 
anisotropic  orientations.  Lead  orientation  and  actual 
width  distributions  are  then  assumed  known,  and 
expected  error  in  lead  width  is  determined  for  a  variety  of 
situations. 

Lead  width  distributions  have  been  described  by  power 
laws  [WodAoms,  1981;  Steffen ,  1987]  as  have  floe  sizes 
[Roth rock  and  Thorndike,  1984].  The  negative  exponential 
distribution  has  also  been  used  [Dirifcins  et  al.,  1986]: 

fyix)  =  e  -%rK  /  X 

with  mean  lead  width  \  and  variance  X2.  The  exponential 
model  implies  that  there  are  a  finite  number  of  small 
leads,  and  that  the  field  is  characterized  by  a  length  scale 
X.  In  fact,  the  lead  width  distribution  may  be  scale- liee, 
in  which  case  a  power  law  would  be  appropriate.  There  is, 
of  course,  a  lower  limit  imposed  by  the  resolution  of  the 
measuring  instrument,  and  for  this  reason  as  well  as  for 
clarity  of  illustrating  expected  values,  we  use  the  negative 
exponential  model. 

Lead  orientations  may  be  random  or  may  have  a 
preferred  orientation.  A  Gaussian  model  is  used  here  for 
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TABLE  1.  Expected  Error  in  Lead  Width*  (in  Meter*)  Under  a  Variety  of  Assumed  Distribution* 

and  Mean  Values 


Case 

fx 

fx 

fe 

♦r 

EOi-X) 

VariX-Jf) 

1 

7 

Sonar 

Uniform 

— 

66.7 

1466.1 

2 

X=20m 

7 

Gaussian^ 

2.36  (135°) 

3.7 

3.1 

3 

X=20m 

7 

Gaussian^ 

0.52  (30°) 

32.5 

165.9 

4 

L=40m 

7 

Gaussian^ 

2.36  (135°) 

4.8 

6.4 

6 

X=40m 

7 

Gaussian^ 

0.52  (30°) 

36.0 

345.3 

6 

Xs20m 

7 

Uniform 

— 

43.2 

653.1 

7 

X=40m 

7 

Uniform 

— 

64.2 

1391.2 

Question  mark*  refer  to  the  unknown  distribution. 

Width  distribution  model  is  negative  exponential. 

*  Parameters  of  the  Gaussian  model  are  \x=kJ4,t  (45°)  and  o=0.3T  (17°). 


preferred  orientations.  It  is  recognized,  however,  that  the 
actual  shape  of  the  distribution  may  be  bimodal,  where 
large  leads  with  one  orientation  are  intersected  by  smaller 
leads  at  another.  Intersection  angles  of  approximately  28° 
have  been  observed  elsewhere  [. Marko  and  Thomson, 
1977].  This  situation  is  not  obvious  in  Figure  4,  although 
the  distribution  is  not  strictly  Gaussian  either. 

Table  1  lists  the  expected  error  for  a  variety  of  condi¬ 
tions,  where  error  is  defined  by  the  difference  between  the 
actual  and  measured  lead  widths.  Case  1  considers  the 
situation  where  the  apparent  lead  width  distribution  is 
known.  The  apparent  lead  widths  are  based  on  submarine 
sonar  data  recorded  by  the  USS  QUEENFISH  in  August 
of  1970  in  the  central  Canada  Basin  [McLaren,  1989].  Ice 
draft  data  were  measured  by  an  upward-beamed  fatho¬ 
meter  with  a  footprint  diameter  of  approximately  2  7  m 
and  a  vertical  accuracy  of  ±10  cm.  Sequences  of  continu¬ 
ous  points  with  drafts  £30  cm  constitute  leads,  an  example 
of  which  is  given  in  Figure  7.  The  apparent  width  distri¬ 
bution  was  determined  for  a  100  km  section  and  is  shown 
in  Table  2.  Given  a  mean  apparent  lead  width  of  60.6  m, 


Ice  Draft  from  Sonar 
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Fig.  7.  Submarine  sonar  ice  draft  date  for  s  2.5  km  section 
within  the  Canada  Basin  north  of  Alaska.  Lead*  are  defined  as 
continuous  sequences  of  points  with  drafts  no  greater  than  0.3  m 
(dashed  fine};  ax  lead*  occur  in  this  section. 


the  expected  value  of  the  error  is  65.7  m  with  a  variance 
of  1466.1  For  cases  2-7  in  Table  1  the  actual  width 
density  function  is  assumed  known,  and  the  apparent  lead 
width  distribution  is  estimated.  For  cases  2  and  3,  the 
crossing  angle  distributions  are  shown  in  Figure  8,  and 
the  error  distributions  in  Figure  9.  In  the  preferred 
orientation  cases  (2-5)  the  error  means  and  variances  are 
clearly  dependent  upon  transect  orientation. 

6.  Conclusions 

A  methodology  has  been  presented  for  determining 
width  distributions  of  linear  features  from  measurements 
along  a  transect  through  a  network  of  such  features.  Both 
isotropic  and  anisotropic  orientations  of  the  linear  features 
have  been  considered.  In  the  anisotropic  case,  the  orienta¬ 
tion  distribution  of  the  lines  must  be  known.  In  both 
cases  if  the  distribution  of  actual  widths  and  the  orienta¬ 
tions  are  viewed  as  independent  random  variables  then 
the  actual  width  distribution  can  be  determined  from  the 


TABLE  2.  Lead  Width*  (Bin  Midpoint)  and  Number  of 
Leads  per  Bin,  Measured  by  Submarine  Sonar 


Width 

N 

Relative 

Frequency 

20.00 

134 

0.441 

40.00 

63 

0.207 

60.00 

18 

0.059 

80.00 

25 

0.082 

100.00 

22 

0.072 

120.00 

12 

0.039 

140.00 

4 

0.013 

160.00 

4 

0.013 

180.00 

5 

0.016 

200.00 

5 

0.016 

220.00 

3 

0.010 

240.00 

2 

0.007 

260  00 

0 

0.000 

280.00 

0 

0.000 

300.00 

3 

0.010 

320.00 

0 

0.000 

340.00 

1 

0.003 

360.00 

3 

0.010 

18.422 
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Transect  Orientations  2  36  (solid).  0  52  (dashed) 


Crossing  Angle  (radians) 

Pig.  8.  The  distribution  functions  of  crossing  angles  F^.  for  cases 
2  (solid)  and  3  (dashed)  in  Table  1. 


apparent  widths  and  vice  versa.  Furthermore,  if  the 
actual  widths  and  the  orientations  vary  jointly  then  the 
apparent  width  can  be  determined  from  the  joint  distribu¬ 
tion.  The  application  presented  was  to  measurements  of 
sea  ice  leads  made  by  submarine  sonar.  The  width 
distributions  measured  from  sonar  illustrated  that  the 
general  shape  is  similar  to  those  derived  from  satellite 
imagery,  but  that  the  errors  in  widths  can  be  significant. 
Unfortunately,  it  is  not  possible  to  determine  the  actual 
error  in  lead  widths  derived  from  much  of  the  archived 
sonar  data.  However,  if  lead  "climatologies"  can  be 
compiled  for  various  locations  and  seasons,  or  if  significant 
relationships  between  lead  orientation  and  geostrophic 
winds  can  be  developed,  then  at  least  we  can  determine 
the  probable  error. 

Of  course,  the  potential  error  is  not  an  issue  if  adequate 
two-dimensional  data  are  available.  For  example,  side- 
scan  sonar  may  permit  a  more  accurate  retrieval  of  lead 
and  keel  statistics  [ Wadhama ,  1988].  Sonar  data  with 

Error  Distributions 

Transect  Orientations  2  36  (solid)  0  52  (dashed) 


Error  (m) 

Pig.  9.  The  distribution  functions  of  the  lead  width  errors,  Fy  jf, 
for  esses  2  (solid)  and  3  (dashed)  in  Table  1. 


concurrent  overhead  imagery  from  aircraft  or  satellite  is 
also  a  potentially  valuable  source  of  information.  How¬ 
ever,  lead  width  distributions  derived  from  satellite  data 
may  not  be  accurate  since  very  small  leads  are  not 
resolved.  This  is  particularly  true  for  medium  resolution 
data  such  as  that  from  the  Defense  Meteorological  Satel¬ 
lite  Program  (DMSP)  or  the  Advanced  Very  High  Resolu¬ 
tion  Radiometer  (AVHRR)  on-board  the  TIROS-N  satel¬ 
lites.  One  solution  might  be  to  retrieve  the  orientation 
information  from  satellite  data  where  small  lead  widths 
cannot  be  resolved,  and  the  width  information  from 
submarine  sonar  transects. 

Other  applications  of  this  procedure  are  possible.  For 
example,  laser  profilometer  transects  are  analogous  to 
sonar  transects,  and  the  methods  outlined  above  could  be 
used  for  lead  and  ridge  spacing  distributions  and  their 
associated  errors.  As  in  the  illustration  with  Landsat 
data,  transect  sampling  of  satellite  imagery  is  a  natural 
application.  Similarly,  heat  flux  through  leads  is  in  part 
a  ftmction  of  fetch,  and  fetch  is  a  function  of  the  actual 
lead  width  and  the  crossing  angle  of  the  wind.  If  the  wind 
direction  is  constant  as  it  travels  across  the  network  of 
leads,  then  the  distribution  of  fetches  can  be  determined 
from  the  distribution  of  actual  lead  widths.  Finally,  it 
may  be  possible  to  estimate  open  water  fraction  over  a 
large  area  from  the  apparent  lead  width  and  spacing 
distributions  measured  along  a  transect.  This  research  is 
currently  in  progress,  with  results  to  be  presented  subse¬ 
quently. 

Appendix 

Derivation  of  fa  From  /j  X 

Equation  (6)  gives  an  expression  for  However, 

we  are  interested  in  fa  rather  than  /|og  80  we  use  the 
identity 

f\og  yflog  y)  -  y  fY<y~>  (A1) 

which  is  proven  by 

*iog  yfl°8  y)  =  Wog  Y  £  log  y] 

=  P[elogTs  elogy] 

=  P[Y  <  y] 

=  Fyiy) 

The  derivative  of  both  sides  with  respect  to  y  produces 
(Al). 

Using  (Al)  and  (6),  an  expression  for  the  actual  lead 
width  distribution  can  be  obtained: 


flog  X<l°g  *>  =  <P 


yl/iog  yl 

•pt/iog  yl 


(log  l) 


fx(x)  =  1  <p'x 


(A2) 


Matrix  Formulation  of  fa 

Expressing  equation  (5)  in  terms  of  matrices  provides  a 
way  to  solve  for  fa.  The  functions  can  be  discretized  as 
arrays  and  the  integral  approximated  as  a  sum  as  shown 
in  (9),  which  is  repeated  here: 

Fxii)  '  A  £  FAij,i)fx(i),  (9') 

i=l 

»j€[l  An,i-iA,  y-j& 
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where  F^JJ.i)  =  FA<lcoa'Hi  fj)],  N  is  the  number  of  discrete 
observations  and  A  is  the  increment  between  observations. 
Suppose  we  now  form  an  N  x  N  matrix,  M,  as  follows: 

M(J,i)  =  A  FA4J,i)  ,  1  £i<j<N 
M{j,i)  =  0  ,  l<j<iiN 

where  the  first  expression  refers  to  elements  on  or  below 
the  diagonal.  It  follows  that 

j  N 

A  £  FA'(j,i)  fx(i)  =  £  M(j,i)  fx(i)  (A3) 
i*l  «=1 

or 

N 

Fx<J)  =  J2  fx(i ) 

»=1 

This  is  equivalent  to  the  matrix  equation 

Fx  -  M  fx  (A5) 

From  this  it  can  be  seen  that  the  problem  is  solved  if  M  is 
invertible.  However,  since  M(jJ)  *  F^/jJ)  *  0  for  all  j 
(from  the  definition  of  FA>),  the  matrix  M  has  zeros  on  as 
well  as  above  the  diagonal,  therefore  the  first  row  and  last 
column  contain  all  zeros,  so  that  det(Af)  =  0,  or  M  is  not 
invertible. 

This  difficulty  is  overcome  by  defining  a  new  N- 1  xN-1 
matrix  to  be  the  submatrix  formed  by  deleting  the  first 
row  and  last  column  from  M.  This  gives 


V 

/  \ 

f2 

h 

MM  ‘I 

Fa 

fN-h 

~  J#2 

Fn 

K  n) 

where  f%  -  f?-  ’  frf  and  Fx  m  f2’-’  fn^-  Since 

Fj  s  0,  it  is  not  needed  to  find  Furthermore,  if  N 

and  A  are  chosen  properly,  then  fN  =  0  because  f\ s  a  pdf. 
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Arctic  Ice  Surface  Temperature  Retrieval 
from  AVHRR  Thermal  Channels 
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The  relationship  between  AVHRR  thermal  radiances  and  the  surface  (skin)  temperature  of 
Arctic  snow-covered  sea  ice  is  examined  through  forward  calculations  of  the  radiative  transfer 
equation,  providing  an  ice/snow  surface  temperature  retrieval  algorithm  for  the  central  Arctic 
Basin.  Temperature  and  humidity  profiles  with  cloud  observations  collected  on  an  ice  island 
during  1986-1987  are  used.  Coefficients  that  correct  for  atmospheric  attenuation  are  given  for 
three  Arctic  clear  sky  ’seasons',  as  defined  through  statistical  analysis  of  the  daily  profiles,  for  the 
NOAA  7,  9,  and  11  satellites.  Modeled  directional  snow  emissivitiee,  different  in  the  two  split- 
window  (11  and  12  urn)  channels,  are  used.  While  the  sensor  scan  angle  is  included  explicitly  in 
the  correction  equation,  its  effect  in  the  dry  Arctic  atmosphere  is  small,  generally  leas  than  0.1  K. 
Using  the  split-window  channels  and  scan  angle,  the  rms  error  in  the  estimated  ice  surface 
temperature  (1ST)  is  leas  than  0.1  K  in  all  seasons.  Inclusion  of  chsnnel  3  (3.7  pm)  during  the 
winter  decreases  the  rms  error  by  lass  than  0.003  K.  The  seasonal  dependence  of  the  coefficients 
is  important,  with  errors  in  the  range  of  0.1  -  0.6  K  when  coefficients  from  one  season  are  used 
with  data  from  another.  Similarly,  mixing  coefficients  and  data  from  different  satellites  results 
in  average  errors  from  0.1  to  1.0  K.  Overall,  employing  the  1ST  coefficients  results  in  increased 
accuracy  of  up  to  0.8  K  over  SST  coefficients  developed  for  the  North  Atlantic  and  the  Greenland 
Sea  areas. 


1.  Introduction 

Satellite  data  for  the  estimation  of  radiative  and 
turbulent  h«at  fluxes  is  becoming  an  increasingly  impor¬ 
tant  tool  in  large  scale  studies  of  climate.  One  parameter 
needed  in  the  estimation  of  these  fluxes  is  surface  temper¬ 
ature.  Sea  and  land  surface  temperature  (SST  and  LST) 
retrieval  algorithms  have  been  developed  by  using  the 
thermal  infrared  window  portion  of  the  spectrum,  with  the 
degree  of  success  dependent  primarily  upon  the  variability 
of  the  surface  and  atmospheric  characteristics.  The 
general  approach  to  estimating  surface  temperature  is  to 
relate  satellite  observations  to  surface  temperature 
observations  with  a  regression  model.  Lacking  sufficient 
observations,  however,  satellite  radiances  or  brightness 
temperatures  can  be  modeled  by  application  of  the  radia¬ 
tive  transfer  equation.  This  approach  is  commonly  used 
for  SST  retrieval  [cf.  Minnett,  1990;  Llewellyn- Jones  et  al. , 
1984;  Barton,  1986],  A  more  complete  review  of  SST 
algorithms  is  given  by  Me'  -in  et  al.  11985].  For  SST 
estimated  using  two  "split- window"  infrared  channels  (e.g., 
approximately  1  pm  wide  centered  at  10.8  and  12.0  pm) 
an  absolute  accuracy  of  0.6-1  K  (rms  error)  has  been 
obtained  \JJewellyn-Jones  et  al.,  1984;  McClain  et  al., 
1985].  Land  surface  temperature  estimation  is  generally 
less  accurate  due  to  the  larger  variability  of  surface 
conditions,  where  errors  of  *2-3  K  are  common  [Price, 
1983]. 

To  our  knowledge,  little  effort  has  been  directed  to  the 
retrieval  of  the  sea  ice  surface  temperature  (1ST)  in  the 
Arctic,  an  area  where  the  first  effects  of  a  changing 
climate  are  expected  to  be  seen.  The  reason  is  not  one  of 
methodology,  but  rather  our  limited  knowledge  of  atmo¬ 
spheric  temperature,  humidity,  and  aerosol  profiles,  the 
microphysical  properties  of  polar  clouds,  and  the  spectral 


characteristic*  of  the  surface  types  found  there.  In  this 
paper  we  present  a  means  to  correct  for  the  atmospheric 
attenuation  of  satellite-measured  clear  sky  brightness 
temperatures  used  in  the  retrieval  of  snow-covered  ice 
surface  temperature  from  the  split-window  thermal 
channels  of  the  advanced  very  high  resolution  radiometer 
(AVHRR)  sensors  on  board  three  of  the  NOAA  series 
satellites.  These  corrections  are  specified  for  three 
different  "seasons"  <»nd  as  a  function  of  satellite  viewing 
angle  and  are  expected  to  be  applicable  to  the  perennial 
ice  pack  in  the  central  Arctic  Basin  (Figure  1).  We  dA  not 
develop  a  completely  new  methodology;  instead  we  rr.  dify 
a  standard  procedure  "or  use  with  Arctic  data.  In  this 
paper  it  is  assumed  that  a  valid  cloud-clearing  algorithm 
exists  and  that  only  clear  sky  radiances  are  being  exam¬ 
ined.  The  cloud  clearing  problem  in  polar  satellite  data  is 
not  trivial,  however.  For  a  review  of  polar  cloud  detection 
algorithms,  see  Key  and  Barry  [1989]  and  Sakellariou  et 
al.  [1992]. 

2.  Data 

Analyses  are  based  on  atmospheric  temperature  and 
humidity  profiles  collected  by  rawinsonde  from  a  Soviet  ice 
island  (NP-26),  located  at  approximately  85°N  170°W 
during  1983-1987  (Figure  1).  Generally,  two  observations 
per  day  were  collected  covering  a  vertical  range  of  0  -  25 
km.  Profiles  that  have  at  least  10  levels  are  retained  in 
the  analysis.  Observations  include  temperature,  dew 
point  depression,  wind  speed,  and  wind  direction.  Ice 
crystal  precipitation  is  not  reported.  For  the  years  1986- 
1987,  surface-based  cloud  observations  are  also  available. 
These  observations  include  low,  middle,  and  high  cloud 
types,  height  of  the  cloud  base,  and  cloud  fraction. 

Only  clear  sky  profiles  are  used  in  this  study  (1986- 
1987),  since  the  satellite  thermal  radiances  under  cloudy 
conditions  will  measure  cloud  top  temperature  and  a 
significant  amount  of  cloud  cover  will  affect  the  lower 
tropospheric  temperature  structure.  Clear  sky  is  defined 
to  be  no  more  than  25%  cloud  cover  by  the  surface  obser¬ 
vations.  Clear  sky  "seasons"  that  differ  in  their  vertical 
temperature  and  humidity  structures  are  then  defined. 
The  seasons  are  determined  objectively  with  a  squared 
Euclidean  distance  clustering  algorithm;  the  variables  are 
temperature  and  humidity  at  each  level.  To  reduce  the 
degree  of  statistical  dependence  between  levels,  only  one 
measurement  per  kilometer  was  used.  The  resulting 
seasons  are  winter  (October  through  March),  summer 
(June  through  August),  and  transition  (April,  May,  and 
September).  This  analysis  was  also  performed  with  three 
other  methods:  subjectively,  clustering  principal  compo¬ 
nent  scores,  and  using  a  simple  correlation  method.  The 
only  difference  between  them  was  in  the  placement  of 
October.  By  temperature  alone  it  is  grouped  vith  the 
transition  months.  When  humidity  profiles  are  also 
included  in  the  analysis,  however,  it  is  more  similar  to  the 
other  winter  months. 

The  resulting  mean  seasonal  temperature  profiles  for 
clear,  cloudy  (greater  than  75%  cloud  cover),  and  mixed 
conditions  are  shown  in  Figure  2.  Differences  reflect  not 


only  the  near-surface  temperatures  but  also  the  vertical 
temperature  structure,  which  can  also  be  seen  in  the 
monthly  inversion  frequencies  reported  by  Serreze  et  al. 
[1991].  The  effect  of  cloud  on  the  surface  radiation 
balance  is  evident  in  the  clear  and  cloudy  profiles  where 
temperatures  in  the  lower  troposphere  are  higher  under 
cloud  cover  in  the  winter  and  transition  but  lower  under 
summer  cloud  cover.  In  addition,  the  summer  surface 
inversion  disappears  under  cloud  cover.  For  comparison, 
the  lowest  30  km  of  the  standard  subarctic  summer  and 
winter  temperature  profiles  from  the  LOWTRAN  7 
radiative  transfer  model  data  base  [Kneizys  et  al.,  1988] 
are  shown  in  Figure  3.  The  most  obvious  differences 
between  the  subarctic  and  Arctic  profiles  are  the  surface 
temperature  and  low-level  inversion  structure.  In  the 
Arctic  profiles,  standard  deviation  at  all  levels  was 
computed  for  each  parameter.  The  winter  standard 
deviations  of  tropospheric  clear  sky  temperature  and  dew 
point  depression  were  computed  for  each  layer,  between  0 
and  10  km.  The  maximum  values  occurred  near  the 
surface  and  were  6.8  K  and  2.7  K,  respectively.  In 
summer  the  corresponding  values  are  4.5  K  and  6.7  K, 
and  in  transition  they  are  7.9  K  and  4.2  K.  Interpretation 
of  dew  point  depression  variability  requires  care  since  the 
mean  water  vapor  content  is  very  low,  generally  less  than 
0.6  gm  m"3. 

The  AVHRRs  on  board  the  NOAA  7, 9,  and  11  satellites 
are  of  interest  in  this  study.  Of  the  five  AVHRR  channels 
the  three  thermal  channels  (channel  numbers  3,  4,  and  5 
centered  at  approximately  3.7, 11,  and  12  jim)  are  simulat¬ 
ed,  and  actual  data  are  used  for  validation.  NOAAs  8  and 
10  are  not  used  because  they  lack  channel  5.  First-order 
calibration  was  performed  following  the  methods  described 
by  NOAA  [1991]  and  Lauritsen  et  al.  [1979],  Additional 
corrections  were  applied  to  the  data  to  account  for  the 
nonlinear  response  of  the  thermal  channels  [  Weinreb  et  al. , 
1990].  Only  NOAA  7  and  11  data  were  used  in  the 
validation.  The  AVHRR  scan  angle  ranges  from  0°  to 
app.  o*.  nately  55°.  Both  local  area  coverage  (LAC;  1.1-km 
field  of  view  at  nadir)  and  global  area  coverage  (GAC; 
approximately  4  km)  data  are  used  in  validation. 

3.  Retrieval  Methodology 

For  the  retrieval  of  SST  a  multichannel  algorithm  that 
corrects  for  atmospheric  attenuation  of  upwelling  radiation 
primarily  due  to  water  vapor  absorption  is  commonly 
employed  [e.g.,  Barton  et  al.,  1989]: 

T„  *  <K 6)  +  £  6/6)7;  <D 

«-3 

where  a(6)  and  6/6)  are  satellite  zenith  angle-dependent 
coefficients  and  Tt  are  the  satellite-measured  brightness 
temperatures  in  the  three  AVHRR  thermal  channels.  The 
coefficients  are  determined  through  a  least  squares 
regression  procedure,  where  surface  temperatures  are 
regressed  against  modeled  brightness  temperatures.  It  is 
also  possible  to  use  differences  between  or  ratios  of  two 
channels.  Such  an  approach  was  taken  by  Schlueeeel  and 


G-aaal  [1990]  for  SST  retrieval  at  high  latitudes. 

Alternatives  to  computing  a  different  set  of  coefficients 
for  each  scan  angle  increment,  as  shown  in  (1),  were 
sought.  Equations  that  include  explicitly  the  scan  angle 
have  been  presented  [cf.  McClain  et  al.,  1986].  Here  we 
use  the  equation 

=  a  *  +  cT5  *  d[(T4-T6Uecel  (2) 

for  1ST  retrieval,  where  brightness  temperatures  are  in 
Kelvin.  The  coefficients  are  given  in  the  next  section.  The 
overall  equation  (2)  produced  the  smallest  standard  error 
of  all  combinations  of  channels,  channel  differences,  and 
scan  angle  functions  tested,  e.g.,  sec  9,  [(1-sec  OXT^-Tg)], 
and  similar  variations.  Although  channel  3  could  be 
included  in  (2),  its  use  would  be  limited  to  winter  analy¬ 
ses,  since  it  measures  reflected  solar  radiation  as  well  as 
thermal  emissions.  The  usefulness  of  this  channel  in  1ST 
retrieval  is  discussed  in  the  next  section. 

To  simulate  radiances  in  the  AVHRR  thermal  channels, 
the  daily  temperature  and  humidity  profiles  in  each 
season  are  used  with  the  LOWTRAN  7  (hereafter 
LOWTRAN)  radiative  transfer  model  [ Kneizys  et  al. ,  1988]. 
Earlier  versions  of  LOWTRAN  have  been  used  in  the 
retrieval  of  SST  [e.g.,  Barton  et  al.,  1989],  LOWTRAN 
calculates  atmospheric  transmittance/radiance  for  wave 
numbers  ranging  from  0  to  60,000  cm'1  (0.2  pm  to  infinity) 
at  a  resolution  of  20  cm'1  (for  gaseous  absorption)  and 
includes  calculations  for  multiple  scattered  radiation.  The 
code  may  be  initialized  for  standard  or  user-defined 
atmospheres,  several  cloud  models,  aerosol  models,  and 
specified  solar  and  view  geometTy.  Radiances  are  calculat¬ 
ed  at  5  cm'1  intervals  (interpolated  by  LOWTRAN), 
equivalent  to  0.06  pm  at  11  pm.  In  this  study,  radiances 
are  modeled  for  sensor  scan  angles  from  0°  to  60°  in  10° 
increments.  Atmospheric  chemical  composition  and  back¬ 
ground  tropospheric  and  stratospheric  aerosols  for  the 
subarctic  winter  and  summer  models  are  used,  since  no 
such  information  is  available  from  the  ice  islands.  The 
optical  properties  of  Arctic  haze  have  not  been  extensively 
measured;  model  calculations  [Blanchet  and.  List,  1983] 
show  that  the  volume  extinction  coefficient  of  Arctic  haze 
is  generally  of  the  same  order  of  magnitude  as  that  of  the 
tropospheric  aerosols  [Taay  et  al. ,  1989].  Therefore  the  use 
of  tropospheric  background  aerosols  is  appropriate.  The 
appropriate  sensor  response  function  (Figure  4)  is  applied 
to  the  calculated  radiances,  and  radiances  are  then 
converted  to  brightness  temperatures  by  inverting  the 
Planck  function  at  the  central  wavelength  of  the  channel 
appropriate  for  the  temperature  range  230  -  270  K  [NOAA, 
1991].  The  appendix  provides  additional  detail  on  simu¬ 
lating  satellite  radiances. 

The  ice  surface  is  assumed  to  be  snow-covered  year 
round.  Directional  surface  emissivities  for  snow  are 
modeled  following  the  procedure  of  Dozier  and  Warren 
[1982].  Briefly,  the  single  scattering  albedo  and  asymme¬ 
try  factor  in  the  scattering  phase  function  are  calculated 
from  the  Mie  equations,  and  the  directional,  wavelength- 
dependent  emissivities  are  derived  from  the  delta- 
Eddington  approximation  to  the  equation  of  radiative 


transfer.  The  directional  emissivities  are  then  integrated 
with  the  response  function  for  channel  i: 

C*  t(X,e)<fc<X)  dx 


where  r(X,0)  is  the  emissivity  in  direction  0  at  wavelength 
X  and  p  is  the  sensor  response  function  which  is  0  outside 
of  [X^].  These  emissivities  are  given  in  Table  1  for 
NOAA  7.  The  difference  between  the  integrated  emissivit¬ 
ies  for  the  three  satellites  examined  here  is  approximately 
0.0001,  so  we  only  use  those  modeled  for  NOAA  7.  At  the 
wavelengths  of  AVHRR  channels  4  and  6  the  emissivity  is 
essentially  insensitive  to  snow  graLi  size  as  well  as  the 
amount  of  liquid  water,  up  to  20%  of  the  total  particle 
volume.  Additionally,  these  emissivities  do  not  change 
significantly  over  the  range  of  temperatures  encountered 
and  therefore  are  applicable  to  the  snow  types  encoun¬ 
tered  year  round.  It  is  possible,  however,  that  melt  ponds 
will  contaminate  the  field  of  view  during  the  summer.  We 
do  not  adjust  emissivities  to  account  for  this  phenomenon. 

The  use  of  the  rawinsonde  profiles  in  modeling  the 
surface  temperature  requires  an  additional  step,  since  the 
first  measurement  in  each  profile  is  the  shelter  tempera¬ 
ture,  not  the  surface  temperature.  (In  the  discussions  that 
follow,  the  terms  "skin"  and  "surface"  temperature  are 
used  interchangeably.)  The  shelter-surface  temperature 
difference  can  be  significant:  more  than  10  depending  on 
the  region  and  time  of  year  [cf.  Stowe  et  ai.,  1988;  Rossow 
et  al.,  1989].  Therefore  the  (unknown)  surface  tempera¬ 
ture  for  each  profile  is  assigned  a  series  of  values  repre¬ 
senting  the  range  of  possible  surface  temperatures  for  the 
observed  conditions  during  the  month  to  which  the  profile 
belongs.  An  energy  balance  model  is  used  to  determine 
these  surface  temperatures,  based  on  the  observed  range 
of  shelter  temperatures  and  wind  speeds  (the  mean  ±1 
standard  deviation)  in  the  ice  island  data  for  each  month. 
After  Maykut  [1982],  the  energy  balance  equation  is 

(l-a)Fr-/jce*F£*eoT*e+F,*Fe+Fe  =  0 

where  o  is  the  albedo,  e  is  the  longwave  emissivity,  o  is 
the  Stefan-Boltzmann  constant  (in  W  m‘2  K"*),  /ic<  is  the 
amount  of  shortwave  energy  that  penetrates  the  ice  and 
does  not  directly  heat  the  surface,  Fr  and  Fj  are  the 
downwelling  shortwave  and  longwave  radiation  respective¬ 
ly,  F$  and  F(  are  the  sensible  and  latent  heat  fluxes, 
respectively,  and  Fc  is  the  conductive  heat  flux.  A  flux 
toward  the  surface  is  positive.  The  sensible  and  latent 
heat  fluxes  are  dependent  in  part  upon  the  wind  speed,  air 
temperature,  7jc#,  and  bulk  transfer  coefficients.  Radia¬ 
tive  fluxes  are  modeled  based  on  the  mean  monthly 
temperature  and  humidity  ice  island  profiles.  Ice  thick¬ 
nesses  are  taken  from  Maykut  [1982].  Three  new  profiles 
are  created  for  each  original  clear  sky  profile,  where  the 
surface  temperatures  are  the  minimum,  mean,  and 
maximum  equilibrium  temperatures  estimated  by  the 
model.  This  results  in  seasonal  sample  sizes  of  750,  123, 


and  225  profiles  for  the  winter,  summer,  and  transition 
seasons,  respectively. 

4.  Discussion 

Coefficients  for  the  estimation  of  1ST  are  given  in 
Tables  2-4  for  NOAA  7,  9,  and  11  satellites  and  are  used 
with  (2).  Coefficients  are  given  to  five  significant  digits. 
(The  rms  error  for  using  three  significant  digits  rather 
than  six  in  winter,  for  example,  is  0.07  K,  0.003  K  for  four 
digits,  and  0.0007  for  five  digits.)  In  all  cases  the  coeffi¬ 
cient  of  determination  (R2)  is  at  least  0.98.  Also  shown  is 
the  root  mean  square  (rms)  error  for  the  difference 
between  the  actual  (energy  balance  model)  surface  temper¬ 
ature  and  the  satellite-derived  surface  temperature. 
These  errors  are  small — always  less  vtaan  0.1  K  and  some¬ 
times  half  that  value — as  a  result  of  (1)  the  low  water 
vapor  content  of  the  atmosphere  and  therefore  little 
atmospheric  attenuation,  (2)  the  radiative  transfer  model¬ 
ing  approach,  which  does  not  incorporate  satellite  or 
ground-based  instrument  noise,  and  (3)  natural  variability 
in  the  polar  atmosphere  that  was  not  captured  in  the  ice 
island  profiles.  While  the  Arctic  atmosphere  is  very  dry, 
however,  an  atmospheric  correction  as  in  equation  (2)  is 
important.  The  uncorrected  channel  4  and  5  brightness 
temperatures  can  be  significantly  different  from  the 
surface  temperature;  e.g.,  the  mean  T— T4  difference  and 
its  standard  deviation  (summer)  are  0.66  and  0.43  K,  0.96 
and  0.61  K  for  Ticg-T6,  -0.21  and  0.32  K  for  the  winter 
channel  4  difference,  and  -0.23  and  0.43  K  for  channel  6 
in  winter,  based  on  NOAA  7  simulated  radiances  at  a 
surface  emissivity  of  unity.  Maximum  differences  are  in 
the  range  of  1.7  to  2.3  K  in  summer  and  -2.5  to  -3.5  in 
winter,  for  channels  4  and  6,  respectively. 

Due  to  a  lack  of  Arctic  surface  and  atmospheric  data 
over  the  ice,  it  is  difficult  to  define  the  area  for  which 
these  coefficients  apply.  At  present  there  is  no  evidence 
that  temperature  and  humidity  characteristics  differ 
significantly  over  different  portions  of  the  perennial  ice 
pack,  but  this  cannot  be  confirmed.  It  has  been  found, 
however,  that  synoptic  activity  is  quite  different  in  the 
eastern  Arctic  (e.g.,  Kara  Sea),  affecting  the  strength  and 
frequency  of  low-level  inversions  as  well  as  humidity 
profiles  [Serreze  et  al.,  1992].  We  therefore  consider  these 
coefficients  to  apply  to  pack  ice  in  the  central  Arctic  Basin. 
The  applicability  of  these  coefficients  to  Antarctica  is 
uncertain.  The  Arctic  and  Antarctic  are  both  character¬ 
ized  by  low  temperatures,  surface  inversions,  snow-covered 
surfaces,  and  low  water  vapor  amounts.  However,  the 
temperatures  and  water  vapor  amounts  are  lower  over 
Antarctica,  although  it  is  unknown  how  much  lower  the 
water  vapor  amounts  are,  by  virtue  of  its  higher  elevation 
end  lower  incidence  of  cyclonic  systems.  Given  this,  one 
would  expect  there  to  be  some  difference  between 
Antarctic  surface  temperatures  estimated  using  the 
coefficients  presented  here  and  coefficients  based  on 
Antarctic  temperature  and  humidity  profiles,  probably  on 
the  order  of  0. 1  •  0.3  K  (see  the  discussion  of  SST  coeffi¬ 
cients  below). 

The  utility  of  including  channel  3  in  surface  tempera- 


ture  retrieval  has  been  shown  to  be  useful  under  certain 
conditions.  For  example,  LUioellyn-Jonea  et  al.  [1984] 
found  that  triple-window  simulations  for  the  tropics  were 
significantly  better  than  split  window,  but  not  for  temper¬ 
ate  latitudes.  Barton  [1986]  found  channel  3  useful  in 
both  tropical  and  mid-latitude  (Australia)  locations.  For 
1ST  retrieval,  the  use  of  channel  3  would  be  limited  to 
winter  analyses,  since  it  has  a  reflected  solar  component 
as  well  as  the  thermal  emissions.  However,  this  channel 
is  often  noisy,  especially  during  winter,  when  the  amount 
of  energy  emitted  at  those  wavelengths  is  small.  To  test 
its  usefulness,  winter  is  redefined  as  November  • 
February,  thereby  avoiding  a  significant  solar  component 
in  the  upwelling  radiance.  Including  channel  3  in  (2) 
reduces  the  ms  error  by  no  more  than  0.003  K  or  all 
satellites,  which  we  do  not  consider  a  significant  improve¬ 
ment  in  accuracy. 

4.1.  Validation 

Validating  the  coefficients  is  difficult  due  to  the  lack  of 
clear  sky  skin  temperature  measurements  with  corre¬ 
sponding  satellite  data.  Therefore  we  use  both  measured 
and  inferred  surface  temperatures.  For  example,  AVHRR 
GAC  data  over  the  Barents  Sea  during  July  1984 
(NOAA  7)  were  used,  with  the  surface  temperature 
assumed  to  be  near  273.16  K.  This  is  a  reasonable 
assumption  for  melting  snow  but  may  be  an  overestimate 
for  ice  due  to  its  higher  salinity.  Over  pack  ice  near  the 
North  Pole  the  mean  estimated  1ST  was  273.04  K.  For  a 
sample  of  pixels  near  the  marginal  ice  zone  with  some 
melt  ponds  the  mean  1ST  was  273.27  K.  Similarly, 
estimated  ISTs  in  NOAA  11  AVHRR  data  over  Greenland 
and  Baffin  Bay  for  July  1990  averaged  272.82  K  over  the 
ice  sheet  and  272.9  K  over  sea  ice.  The  ice  sheet  location 
examined  was  the  site  of  a  Swiss  Federal  Institute  of 
Technology  (ETH)  camp,  which  reported  melt  conditions 
(K.  Steffen,  personal  communication,  1991). 

GAC  data  during  January  1984  north  of  Greenland 
were  also  used,  where  the  estimated  ISTs  were  compared 
to  temperatures  measured  by  drifting  buoys  [Colony  and 
MuAom,  1986].  There  has  been  some  discussion,  although 
no  formal  study,  concerning  the  accuracy  of  these  buoy 
temperatures.  During  summer  the  buoy  housing  may 
experience  radiational  heating,  and  during  winter  they 
may  be  insulated  by  drifting  snow.  Comparison  of  esti¬ 
mated  ISTs  for  10  pixels  around  the  location  of  two 
different  buoys  on  January  7  yields  mean  temperature 
differences  of  6  K  for  one  buoy  and  11  K  for  the  other. 
The  estimated  ice  surface  temperatures  were  within  1  K 
of  each  other  for  each  set  of  pixels.  Differences  of  this 
magnitude  were  also  reported  by  Com  iso  [1983,  his  Figure 
3]  for  surface  temperatures  estimated  from  temperature 
humidity  infrared  radiometer  (THIR)  data,  so  that  the  use 
of  the  drifting  buoys  for  validation  does  not  appear  to  be 
useful. 

Last,  surface  temperature  measurements  taken  by  a 
PRT  6  thermal  radiometer  during  C  EAR  EX  in  March  1989 
are  compared  to  NOAA  11  AVHRR  data.  The  PRT  6  was 
flown  on  the  NOAA  P3  aircraft  near  Svalbard  as  part  of 


the  Arctic  Gat  and  Aerosol  Sampling  Program  (AGASP). 
The  altitude  of  the  instrument  varied  between  50  m  and 
4  km  for  the  "clear"  sky  area  of  coincident  aircraft/satellite 
data.  While  in  concept  this  data  should  be  useful  for 
validation  of  the  ISTs,  it  is  problematic  due  to  (1)  time 
differences  between  the  aircraft  flight  and  the  satellite 
overpass  (a  few  hours),  (2)  an  assumed  unit  emissivity  in 
the  calculation  of  PRT  5  temperatures,  (3)  geolocation 
errors  for  both  the  aircraft  and  the  satellite  data,  and  (4) 
the  presence  of  aerosols  and/or  ice  crystals  above  the 
aircraft.  The  last  of  these  conditions  greatly  reduced  the 
number  of  areas  usable  for  validation.  The  geolocation 
problem  dictates  that  the  comparison  between  the  PRT  5 
and  AVHRR  data  be  done  over  a  number  of  pixels  rather 
than  a  single  pixel.  Given  these  problems,  the  best 
situation  occurred  when  the  aircraft  was  at  an  altitude  of 
160  m.  The  mean  1ST  for  a  sample  of  four  AVHRR  pixels 
was  258.9  K,  while  the  mean  PRT  5  temperature  (adjusted 
for  an  emissivity  of  0.998)  of  four  consecutive  measure¬ 
ments  1  km  apart  was  259.04  K.  Given  the  difficulties  in 
comparing  the  two  data  sets,  these  results  are  encourag¬ 
ing. 

4.2.  Dependencies  and  Atmospheric  Considerations 

A  potential  problem  in  the  retrieval  of  1ST  is  the 
presence  of  ice  crystal  haze,  or  "diamond  dust."  It  is 
particularly  difficult  to  detect  in  AVHRR  data  because  it 
is  usually  close  to  the  surface  and  exhibits  similar  spectral 
properties.  Values  of  the  visible  optical  depth  for  ice 
crystal  haze  have  been  reported  to  range  from  5  to  21  for 
wintertime  and  from  0.03  to  3  for  springtime  [Curry  et  al. , 
1990].  The  effect  of  varying  optical  depth  of  ice  crystal 
haze  on  the  estimated  1ST  is  illustrated  indirectly  in 
Figure  6,  which  shows  brightness  temperature  differences 
between  AVHRR  channels  4  and  5  as  a  function  of  the 
optical  depth  of  ice  crystal  haze  over  two  ice  surfaces:  6- 
cm-thick  ice  with  a  surface  temperature  of  256.62  K  and 
2  m  ice  with  a  surface  temperature  of  235.38  K.  Mean 
January  atmospheric  conditions  for  the  central  Arctic  are 
used  at  a  satellite  scan  angle  of  30°.  The  top  of  the  ice 
crystal  haze  layer  is  near  the  top  of  the  inversion,  which 
has  a  temperature  of  248  K.  Using  the  coefficients 
developed  here,  estimated  ISTs  over  the  2-m  ice  surface 
for  the  points  indicated  (in  order  of  increasing  optical 
depth)  are  235.36,  237.51, 239.57,  and  246.16  K.  Depend¬ 
ing  on  the  cloud  detection  algorithm  used,  a  temperature 
change  of  more  than  2°-3°  would  probably  signify  cloud,  so 
that  in  this  example,  diamond  dust  with  optical  depths 
greater  than  0.34  do  not  present  a  problem  in  1ST  retriev¬ 
al.  At  smaller  optical  depths,  however,  estimated  ISTs 
could  be  in  error.  Given  the  uncertainty  in  the  frequency 
of  occurrence  and  spatial  extent  of  this  phenomenon,  we 
do  not  attempt  to  adjust  for  it  in  the  coefficients  present¬ 
ed.  Admittedly,  this  can  be  an  important  problem  at 
certain  times  of  the  year. 

The  dependence  of  the  coefficients  on  sensor  scan  angle 
has  been  found  to  be  important  by  other  investigators  [cf. 
Barton,  1985;  Minnett,  1990].  This  is  also  the  case  for  1ST 
retrieval,  although  incorporating  both  channels  4  and  5  in 


(2)  reduces  the  effect  that  increased  path  length  at  large 
scan  angles  has  on  the  surface  temperature  estimation 
when  scan  angle  is  not  taken  into  account  explicitly.  With 
the  coefficients  presented  here  the  estimated  surface 
temperature  from  (2)  is  very  close  to  the  surface  tempera¬ 
ture  at  all  scan  angles,  while  the  estimated  surface 
temperature  using  only  channels  4  and  5  varies  as  a 
function  of  scan  angle.  Even  though  the  differences  are 
small,  not  including  the  scan  angle  explicitly  results  in  an 
increased  rma  error  of  less  than  0.1  K  over  using  (2). 

Coefficients  in  (2)  were  also  computed  assuming 
emissivities  independent  of  wavelength  and  scan  angle.  In 
this  case  a  different  set  of  coefficients  was  produced  at 
emissivities  of  0.96  to  1.0  in  increments  of  0.01.  Using  the 
scan  angle  and  wavelength  dependent  coefficients  (Tables 
2-4)  with  data  based  on  fixed  emissivity  coefficients  results 
in  underestimates  of  0.2  K  for  an  emissivity  of  1.0  and  0.6 
K  for  an  emissivity  of  0.99  during  summer  with  NOAA  7 
simulated  data.  Since  pure  water  has  a  slightly  lower 
emissivity  than  pure  snow  (e.g.,  0.992  at  11.9  pm  and 
nadir),  the  estimated  1ST  of  an  AVHRR  pixel  contaminat¬ 
ed  by  summer  meltponds  would  be  incorrect  by  an  amount 
somewhere  between  these  two  extremes  depending  on  the 
proportion  of  water  within  the  field  of  view. 

The  seasonal  dependence  of  the  coefficients  is  illustrat¬ 
ed  in  Table  6,  where  coefficients  from  each  season  were 
applied  to  data  from  every  other  season.  Results  are 
shown  for  NOAA  9  and  indicate  errors  between  0.1  K  for 
transition  coefficients  with  winter  data  and  0.6  K  when 
summer  coefficients  are  used  with  winter  data.  Similarly, 
the  satellite  dependence  of  the  coefficients  is  shown  in 
Table  6  for  summer  conditions.  On  the  average,  errors 
ranging  from  0.1  to  1.0  K,  depending  on  season,  can  be 
expected  in  applying  coefficients  derived  for  one  satellite 
to  data  from  another,  the  smallest  errors  occurring 
between  NOAA  7  and  9  coefficients  and  data. 

Using  SST  coefficients  developed  for  the  North  Atlantic 
[Llewellyn-  Jones  et  al.,  1984]  and  the  Greenland  Sea  area 
[Minnett,  1990]  to  estimate  1ST  would  result  in  an  under¬ 
estimate  of  up  to  0.7  K,  largest  in  winter  and  at  scan 
angles  of  40°  and  greater.  Not  surprisingly,  the  difference 
is  much  larger  (up  to  a  5.0  K  overestimate)  when  tropical 
coefficients  [Llewellyn- Jones  et  al.,  1984]  are  used.  This 
comparison  was  made  by  using  simulated  radiances  based 
on  the  directional  emissivities,  whereas  the  emissivity  of 
the  sea  surface  in  the  aforementioned  studies  was  comput¬ 
ed  by  using  the  Fresnel  equations  with  the  refractive 
index  of  water  (actual  values  not  given).  These  errors  are 
similar  to  those  reported  by  Minnett  [1986],  where  SSTT 
was  retrieved  from  North  Atlantic  data  using  coefficients 
from  other  regions. 

Sensor  characteristics  and  calibration  also  influence  the 
retrieved  1ST.  AVHRR  GAC  data  are  less  noisy  than  LAC 
data  and  would  result  in  a  smoother  1ST  field.  Incorrect 
sensor  calibration  can  produce  large  errors  in  estimated 
1ST;  errors  in  1ST  near  Greenland  of  up  to  3  K  resulted 
from  not  including  the  nonlinear  response  correction  in 
the  calibration.  Quantization  of  the  signal  (i.e.,  how  much 
of  a  degree  is  represented  by  one  digital  count)  introduces 


additional  error.  Soma  of  these  factors  affect  whether  or 
not  there  is  a  systematic  difference  between  the  estimated 
ISTs  and  the  true  (measured)  1ST — the  system  bias — but 
due  to  the  small  quantity  of  coincident  surface/satellite 
measurements,  the  degree  of  influence  is  impossible  to 
assess.  Therefore  no  empirical  corrections  to  the  forward 
model  are  suggested,  which  would  otherwise  account  for 
sensor  and  other  factors  such  as  radiosonde  accuracy  and 
the  treatment  of  the  stratosphere  above  the  known  profile. 

5.  Summary  and  Conclusions 

The  relationship  between  AVHRR  clear  sky  thermal 
radiances  and  the  surface  (skin)  temperature  of  central 
Arctic  snow-covered  sea  ice  is  examined  through  forward 
calculations  of  the  radiative  transfer  equation.  Tempera¬ 
ture  and  humidity  profiles  and  cloud  data  from  ice  islands 
during  1986-1987  are  used.  Coefficients  that  correct  for 
atmospheric  attenuation  are  given  for  three  Arctic  clear 
sky  "seasons,"  as  defined  through  statistical  analysis  of  the 
daily  profiles,  for  the  NOAA  7,  9,  and  11  satellites. 
Modeled  directional  emissivities,  different  in  the  two  split- 
window  channels,  are  used.  While  the  sensor  scan  angle 
is  included  explicitly  in  the  correction  equation,  its  effect 
in  the  dry  Arctic  atmosphere  is  small,  generally  less  than 
0.1  K.  The  coefficients  presented  are  for  use  with  AVHRR 
channels  4  and  5  (11  and  12  pm)  and  the  sensor  scan 
angle.  With  this  method  the  rms  error  in  the  estimated 
ice  surface  temperature  is  less  than  0.1  K  in  all  seasons. 
Inclusion  of  channel  3  (3.7  pm)  during  the  winter  decreas¬ 
es  the  rms  error  by  less  than  0.003  K 

The  seasonal  dependence  of  the  coefficients  is  impor¬ 
tant,  with  rms  errors  in  the  range  of  0. 1  -  0.6  K  when 
coefficients  are  applied  across  seasons.  Similarly,  using 
coefficients  from  one  satellite  with  data  from  another 
results  in  average  errors  from  0.1  to  1.0  K.  Overall, 
employing  the  1ST  coefficients  results  in  increased  accura¬ 
cy  of  up  to  0.6  K  over  SST  coefficients  developed  for  the 
North  Atlantic  and  the  Greenland  Sea  areas.  While  this 
difference  is  small  in  terms  of  outgoing  longwave  radia¬ 
tion,  it  is  important  for  long-term  climate  monitoring. 
Modeled  ice  crystal  haze  during  January  indicates  that  at 
small  optical  depths  the  haze  may  not  be  detected  as  cloud 
and  could  result  in  1ST  errors  of  up  to  2°. 

Two  important  problems  remain  in  the  retrieval  of  ice 
surface  temperature.  First,  until  a  reliable  method  of 
cloud  clearing  becomes  accepted  by  the  science  community, 
1ST  retrieval  results  will  have  an  additional  level  of 
uncertainty.  Second,  coincident  satellite,  surface,  and 
atmosphere  measurements  over  sea  ice  must  be  taken 
across  the  Arctic  and  in  all  seasons  before  the  bias  of  these 
1ST  coefficients  can  be  accurately  determined. 

Appendix 

Assuming  that  clouds  do  not  contaminate  the  satellite 
scene,  that  the  remaining  atmosphere  is  nonscattering, 
and  that  the  surface  is  nonblack,  then  the  radiative 
transfer  equation  for  the  upward  monochromatic  radiance 
L(X,6)  at  wavelength  X  and  satellite  view  angle  6  can  be 
expressed  as 


L(X,B)  -  ik(B]Bx(Ta)c 

*  jf*  BiSTi&tfe  '% 

*  l(l-ex(«)]«  BxtWx)]*  ^ 


where  is  the  emissivity  of  the  surface  (which  is  as¬ 
sumed  equal  at  all  X  within  a  channel  in  this  study),  Tg  is 
the  surface  temperature  in  Kelvin,  is  the  optical  depth 
of  the  slant  path,  and  B^(T)  is  the  Planck  function  at 
temperature  T.  The  first  term  on  the  right  is  the  contri¬ 
bution  from  the  surface,  the  second  is  from  atmospheric 
emission,  and  the  third  represents  downward  atmospheric 
emission  that  has  been  reflected  upward.  The  surface 
contribution  is  assumed  to  be  the  dominant  one  for  the 
1ST  retrieval  outlined  in  (2),  which  can  be  illustrated  by 
estimating  radiances  in  the  two  split-window  channels 
using  the  clear  sky  Arctic  mean  and  subarctic  standard 
winter  and  summer  profiles  with  identical  surface  temper¬ 
atures.  The  maximum  difference  in  radiances  is  0.05  W 

n  i 

m  ir  ,  indicating  that  the  vertical  temperature  distribu¬ 
tion  of  the  relatively  dry  Arctic  atmosphere  plays  a 
relatively  small  role  in  the  attenuation  of  upwelling 
longwave  radiation. 

To  simulate  the  satellite  radiance,  the  radiances  at 
wavelengths  across  each  channel  must  be  integrated  with 
the  sensor  response  function: 


J L(X,ex.,(X)  dx 
Oi<X)  dX 


where  Li  is  the  channel  i  radiance,  ®j(X)  is  the  channel’s 
response  function,  and  Xj  and  Xj  are  the  lower  and  upper 
limits  of  the  channel,  i.e.,  where  the  response  is  0.  Using 
a  rectangular  response  function  defined  by  the  half¬ 
amplitude  full-width  portion  of  the  actual  channel  re¬ 
sponse  (but  with  100%  response  at  all  wavelengths) 
instead  of  the  full  response  function  results  in  brightness 
temperature  differences  of  the  order  of  0.05  K  in  channel 
4  end  0.6  K  in  channel  5  (January  conditions).  These 
figures  are  valid  for  two  central  wavelengths  used  in  the 
conversion  of  radiance  to  brightness  temperature:  that 
from  the  full  response  function  or  one  based  on  Planck 
radiances  specifically  for  the  half-amplitude  function. 
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TABLE  1.  Angular  Emissivities  of  Snow  in  NOAA  7 
AVHRR  Channels  4  and  5 


Scan  Angle 


Channel  4 
Emissivity 


Channel  5 
Emissivity 


TABLE  2.  Winter,  Transition,  and  Summer  Season  Coefficients  and  rms  Error  Based  on 
AVHRR  Channels  4  and  5  for  NOAA  7 

Season  abed  rms 


Winter 

-3.38668 

6.28508 

-5.27306 

-2.45291 

0.102 

Transition 

-3.77780 

4.73209 

-3.71850 

-1.40115 

0.074 

Summer 

-0.47429 

3.77483 

-2.77389 

-0.56024 

0.067 

Season  coefficients  based  on  equation  (2). 


TABLE  3.  Winter,  Transition,  and  Summer  Season  Coefficients  and  rms  Error  Based  on 
AVHRR  Channels  4  and  5  for  NOAA  9 


Season 

a 

& 

c 

d 

rms 

Winter 

-6.820S9 

7.81491 

-6.79284 

-3.34169 

0.127 

Transition 

-6.06238 

5.64562 

-4.62267 

-1.91927 

0.089 

Summer 

0.49996 

4.12165 

-3.12356 

-0.68087 

0.067 

Season  coefficients  based  on  equation  (2). 
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TABLE  4.  Winter,  Transition,  and  Summer  Season  Coefficients  and  rms  Error  Based  on 
AVHRR  Channels  4  and  5  for  NOAA  11 


Season 

a 

b 

c 

d 

nn b 

Winter 

•5.39436 

5.46800 

•4.45233 

•1.45853 

0.071 

Transition 

-5.35487 

4.47913 

•3.46285 

-0.97128 

0.053 

Summer 

-1.76899 

3.66554 

-2.66249 

-0.39676 

0.053 

Season  coefficients  based  on  equation  (2). 


TABLE  S.  The  rms  Error  in  Applying  Coefficients 
(NOAA  9)  Developed  for  One  Season 
to  Dsts  From  Another 


Data  from 


Coefficients 

Winter 

Summer 

Transition 

Winter 

0 

0.403 

0.128 

Summer 

0.587 

0 

0.342 

Transition 


0.117 


0.219 


0 


TABLE  6.  The  mu  Error  in  Applying  Coefficients 
(Summer)  Developed  for  One  Satellite 
to  Data  From  Another 


Data  from 


Coefficients 

NOAA  7 

NOAA  9 

NOAA  11 

NOAA  7 

0 

0.272 

0.666 

NOAA  9 

0.296 

0 

1.017 

NOAA  11 


0.682 


0.961 


0 


Figures 


Fig.  1.  Average  minimum  (solid)  and  msrimum  (dashed)  sea  ice 
extent  in  the  Arctic  Ocean.  Also  shown  is  the  area  covered  by  the 
NP-26  drifting  ice  island. 

Fig.  2.  Clear,  cloudy,  and  mixed  temperature  profiles  for  the  ice 
islsnd  data  during  winter  (October  -  March),  transition  (April, 
Msy,  September)  and  summer  (June  -  August). 

Fig.  3.  Standard  subarctic  winter  and  summer  temperature 
profiles. 

Fig.  4.  Response  functions  for  AVHRR  channels  4  and  5  on 
NOAA  7,  NOAA  9,  and  NOAA  11. 

Fig.  5.  Brightness  temperature  difference  between  AVHRR 
channels  4  and  5  as  a  function  of  the  optical  depth  of  ice  crystal 
haze  over  two  ice  surfaces:  6-cm-thick  ice  with  a  surface  tempera¬ 
ture  of  263.6  K  and  2-m  ice  with  a  surface  temperature  of  236.4 
K.  Mean  January  atmospheric  conditions  for  the  central  Arctic 
are  used.  The  top  of  the  ice  crystal  haze  layer  is  near  the  top  of 
the  inversion,  which  has  a  temperature  of 246.0  K.  Satellite  scan 
angle  is  30*.  Estimated  ISTs  over  the  2-m  ice  surface  for  the 
points  indicated  (in  order  of  increasing  optical  depth)  are  236.36, 
237.61,  239.67,  and  246.15  K. 


—  -  5  cm  ice 
-  2  m  ice 


I  0.8  1.2 

Ice  Crystal  Optical  Depth 
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ABSTRACT 

The  release  of  heat  from  sea  ice  leads  is  an  important  component  of  the  heat  budget  in  the  Arctic, 
but  their  impact  on  regional  scale  climate  is  difficult  to  assess  without  more  information  on  their 
distribution  in  both  space  and  time.  Remote  sensing  of  leads  using  satellite  data,  specifically  AVHRR 
thermal  and  Landsat  visible  imagery,  is  examined  with  respect  to  one  lead  parameter:  lead  width. 
Atmospheric  effects  are  illustrated  through  the  concept  of  thermal  contrast  transmittance,  where  the 
brightness  temperature  contrast  between  leads  of  various  ice  thicknesses  and  the  surrounding  multiyear 
ice  are  simulated  using  a  radiative  transfer  model.  Calculations  are  made  as  a  function  of  aerosol,  ice 
crystal  precipitation,  and  cirrus  cloud  optical  depths.  For  example,  at  ice  crystal  optical  depths  of  more 
than  0.3  under  mean  January  conditions  in  the  central  Arctic,  the  split-window  brightness  temperature 
differences  between  2  m  and  5  cm  thick  ice  are  similar  to  the  ice  temperature  variability  so  that  there 
would  be  insufficient  contrast  to  distinguish  a  lead  from  the  surrounding  ice.  The  geometrical  aspects  of 
the  sensor  are  also  simulated  so  that  the  effect  of  sensor  field-of-view  on  retrieved  lead  width  statistics 
can  be  assessed.  This  is  done  using  Landsat  data  degraded  to  AVHRR  pixel  sizes. 

INTRODUCTION 

Sea  ice  leads  (linear  openings  in  the  ice  pack)  are  an  important  component  in  local  scale  heat 
exchange  in  the  Arctic,  providing  a  significant  source  of  heat  and  moisture  to  the  atmosphere. 
Unfortunately,  their  impact  on  regional  scale  climate  is  difficult  to  assess  without  more  information  on 
their  frequency  of  occurrence  in  both  space  and  time.  The  use  of  satellite  data,  and  in  particular  the 
thermal  channels  of  the  Advanced  Very  High  Resolution  Radiometer  (AVHRR),  for  the  detection  and 
mapping  of  leads  has  been  suggested  but  its  utility  is  questionable  because  of  its  relatively  low  resolution 
(1.1  km  at  nadir)  and  the  influence  of  the  intervening  atmosphere  on  upwelling  radiation.  It  is  necessary, 
therefore,  to  understand  bow  sensor  characteristics,  atmospheric  properties,  and  surface  conditions 
influence  the  detection  and  interpretation  of  sea  ice  leads  in  AVHRR  and  other  satellite  data.  We  seek 
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to  determine  the  sources  and  magnitudes  of  errors  inherent  in  the  measurements  and  how  lead  statistics 
change  with  the  different  spatial  resolutions  of  existing  and  future  sensors.  Specifically,  we  address  the 
following  issues:  What  is  actually  being  measured  by  these  sensors?  To  what  extent  do  surface, 
atmosphere,  and  sensor  factors  affect  lead  observations  and  image  interpretations?  How  does  lead 
detection  depend  on  ice  thickness  as  a  function  of  sensor  response,  surface  temperatures,  and  atmospheric 
conditions?  Since  the  temperature  contrast  between  open  water  and  ice  provides  a  means  to  map  leads 
using  thermal  imagery,  to  what  degree  does  this  contrast  affect  the  apparent  width  of  a  lead  and  our  ability 
to  detect  it?  How  do  atmospheric,  boundary  layer,  and  solar  zenith  angle  effects  combine  with  surface 
conditions  and  viewing  geometry  to  determine  the  thresholds  of  lead  detection? 

To  address  these  questions  our  approach  includes  both  modeling  and  empirical  studies.  The 
radiative  transfer  modeling  is  done  for  the  purpose  of  simulating  the  satellite  sensor  response  under  a 
variety  of  atmospheric  and  surface  conditions.  The  empirical  studies  include  a  comparison  of  lead 
statistics  determined  in  imagery  of  varying  scales  (e.g.,  AVHRl< ,  Landsat,  and  OLS  data)  and  the  use  of 
lead  statistics  derived  from  low  resolution  imagery  to  estimate  characteristics  of  the  distributions  obtained 
in  the  high  resolution  images.  In  this  paper  we  discuss  both  the  empirical  and  modeling  studies  with 
respect  to  lead  width,  which  is  defined  as  the  distance  across  a  lead  along  a  line  perpendicular  to  the  local 
orientation  of  the  lead. 

2.  THERMAL  CONTRAST  BETWEEN  LEADS  AND  THE  SURROUNDING  ICE 

Atmospheric  and  view  angle  effects  are  modeled  under  a  variety  of  conditions  to  simulate  sensor 
response  in  order  to  gain  a  better  understanding  of  how  leads  might  appear  in  the  satellite  data  under 
varying  surface  and  atmospheric  conditions.  With  this  know  ledge,  it  may  be  possible  to  estimate  the 
smallest  resolvable  lead  under  a  given  set  of  surface/atmospheric/sensor  resolution  and  view  angle 
conditions.  Emphasis  is  on  the  AVHRR  thermal  channels  (3:  centered  at  approximately  3.7  pm,  4:  1 1  pm. 


5:  12  pm). 
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To  simulate  radiances  measured  by  the  sensor,  the  LOWTRAN  7  (hereafter  LOWTRAN) 
atmospheric  transmittance/radiance  model  ( Kneizys  et  a!.,  1988)  is  used.  LOWTRAN  calculates 
atmospheric  transmittance/radiance  for  wave  numbers  ranging  from  0  to  50,000  cm'1  (0.2pm  to  infinity) 
at  a  resolution  of  20  cm'1,  and  includes  calculations  for  multiple  scattered  radiation.  The  code  may  be 
initialized  for  standard  or  user-defined  atmospheres,  several  cloud  models,  aerosol  models,  and  specified 
solar  and  view  geometry.  Analyses  are  based  on  atmospheric  temperature  and  humidity  profiles  collected 
by  radiosonde  from  the  Soviet  ice  islands  RVAB  (approximately  80°N  160°W)  during  1979-83  and  EMIO 
(approximately  87°N  175°E)  during  1983-1987.  Generally  two  observations  per  day  were  collected 
covering  a  vertical  range  of  0  to  25  km.  Profiles  that  had  more  than  10  levels  were  retained  in  the 
analyses.  Pertinent  observations  include  temperature  and  dew  point  depression.  Results  for  mean  January 
conditions  are  presented  here. 

LOWTRAN  calculations  were  made  ever  a  range  of  0.55  micron  optical  depths  for  3  prescribed 
surface  temperatures  used  to  characterize  different  lead  types  (open  or  refrozen)  and  a  fourth  to  represent 
the  "background"  scene  consisting  of  2  meter  thick  ice.  These  temperatures  were  estimated  using  a  surface 
energy  balance  approach  ( Maykut ,  1982)  to  characterize  the  open  leads  at  271  K,  leads  covered  by  5  cm 
thick  ice  (257  K)  and  15  cm  thick  ice  (248  K).  The  thick  ice  has  a  skin  temperature  of  235  K  in  all 
model  calculations.  The  model  was  run  for  satellite  viewing  angles  of  0°  (nadir),  20°  and  50°  and 
accounted  for  differences  in  directional  emissivity  at  each  central  wavelength  of  the  AVHRR  channels. 
Temperature  and  humidity  profiles  were  based  on  the  mean  January  ice  island  soundings  discussed  above. 
All  simulated  radiances  were  converted  to  brightness  temperatures  for  ease  of  analysis  and 
intercomparisons.  Examples  of  the  results  are  shown  in  Figure  1  for  two  commonly  observed  conditions 
in  the  Arctic,  ice  crystal  precipitation  and  haze  in  the  boundary  layer.  Figure  la  is  a  composite  of  results 
at  50°  viewing  angle  showing  the  behavior  of  channel  brightness  temperatures  (numbered  3,  4  and  5  on 
the  plots)  for  simulations  assuming  each  of  the  surface  temperatures  defined  above.  A  hypothetical  cirrus 
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cloud  of  varying  optical  depth  was  inserted  in  the  lowest  1  km  of  the  atmosphere  to  simulate  the  effects 
of  "ice  crystal  precipitation"  which  is  a  common  feature  affecting  the  radiative  balance  of  the  lower 
troposphere  in  the  Arctic.  Regardless  of  the  underlying  surface  type  the  simulated  brightness  temperatures 
rapidly  converge  to  the  248  K  blackbody  radiating  temperature  of  the  "cloud  top"  level  as  optical  depth 
increases.  Note  that  because  all  three  channels  are  centered  in  window  regions  of  the  spectrum,  surface 
temperatures  are  well  simulated  at  zero  optical  depths.  Maximum  brightness  temperature  differences  exist 
between  the  channels  at  around  optical  depth  0.5  to  0.8.  If  satellite  radiances  can  be  measured  accurately, 
these  differences  may  be  used  to  delineate  leads  when  optically  thin  ice  crystal  precipitation  is  present  in 
the  boundary  layer,  but  brightness  temperature  differences  and  contrast  between  surface  types  is  rapidly 
diminished  as  the  optical  depth  increases  beyond  about  1.0,  making  it  impossible  to  detect  the  underlying 
leads.  This  value  is  not  at  all  unexpected  as  optical  depths  of  up  to  20  have  been  reported  for  ice  crystal 
precipitation  ( Curry  et  al.,  1990). 

Figure  lb  shows  the  results  assuming  that  an  aerosol  (haze)  layer  exists  from  the  surface  to  a 
altitude  of  2  km  (the  boundary  layer  aerosol  model  defined  in  LOWTRAN)  over  the  same  range  of  optical 
depths.  In  this  figure  two  sets  of  channel  3  curves  are  shown  to  contrast  the  differences  between  nighttime 
(strictly  thermal)  and  daytime  simulations  where  the  upper,  bold  set  of  curves  relate  to  daytime  conditions. 
Multiple  scattered  solar  radiation  adds  significantly  to  the  upwelling  radiances  as  optical  depths  of  the  haze 
layer  increase  while  the  infrared  channels  do  not  show  any  multiple  scattering  effects.  All  other  input 
parameters  to  LOWTRAN  are  the  same  as  for  Figure  la.  In  the  haze  simulations  it  is  obvious  that  the 
brightness  temperature  differences  between  channel  pairs  are  small  when  only  thermal  radiances  are 
considered  (nighttime)  over  the  entire  range  of  optical  depths  while  large  differences  exist  between  channel 
3,  and  channels  4  or  5  when  solar  radiation  is  accounted  for  during  the  day.  Also,  the  contrast  between 
surface  types  for  each  channel  remain  significant  over  the  range  of  visible  optical  depths  shown, 
particularly  when  only  thermal  radiances  are  considered.  Gross  differences  between  Figures  la  and  lb 
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are  due  to  markedly  different  properties  of  the  haze  and  cirrus  models  defined  for  LOWTRAN. 

The  contrast  between  varying  lead  types  and  the  ice  pack  can  be  expressed  by  the  difference  in 
simulated  channel  brightness  temperature  for  a  hypothetical  lead  temperature,  say  open  water,  5  cm  thick 
ice  or  15  cm  thick  ice  and  the  surrounding  ice  pack  at  any  optical  depth,  normalized  by  the  brightness 
temperature  of  the  ice  pack  (here  taken  to  be  2  m  thick).  Figures  2a  and  2b  show  the  normalized  contrast 
for  channel  4  assuming  the  conditions  illustrated  in  Figures  la  and  lb  respectively.  This  relative  measure 
of  contrast  may  provide  a  means  to  detect  leads  using  AVHRR  data  even  at  night,  though  it  is  apparent 
that  if  optically  thick  ice  crystal  precipitation  exists  that  the  contrast  is  essentially  zero  making  it 
impossible  to  detect  leads  using  such  a  method. 

Also  modeled,  though  not  illustrated  here,  were  channel  brightness  temperatures,  pair  differences  and 
normalized  contrasts  in  surface  types  for  high-level  cirrus  clouds  at  each  viewing  angle.  Again,  brightness 
temperatures  converged  at  the  blackbody  radiating  temperature  at  cloud  top  for  large  optical  depths,  but 
the  convergence  was  less  dramatic  than  in  the  case  of  low  level  ice  crystal  precipitation.  Channel  3  minus 
channel  4  brightness  temperature  differences  at  50°  viewing  angle  were  significant  in  this  case  reaching 
nearly  12  K  over  open  leads  at  a  0.55  micron  optical  depth  of  about  2.6  with  the  difference  being  less  than 
5  K  over  the  pack  ice.  Differences  in  these  signatures  when  incorporated  into  image  analysis  algorithms 
in  conjunction  with  ancillary  information  should  be  very  useful  in  the  detection  of  leads.  Though  the 
normalized  contrast  between  surface  types  degrades  with  increasing  cirrus  optical  depth,  our  simulations 
indicate  that  reasonable  signatures  should  be  detectable  for  a  range  up  to  about  optical  depth  4.0,  that  is, 
for  most  cirrus  conditions  expected  to  be  observed  in  the  Arctic  region.  Combining  brightness  temperature 
differences  with  contrast  will  improve  our  ability  to  detect  leads  significantly.  Developing  an  algorithm 
to  do  so  presents  a  challenge  for  the  future. 
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3.  EMPIRICAL  APPROACHES  AND  RESULTS 

While  there  have  been  studies  of  the  effect  of  sensor  resolution  on  parameter  retrieval,  the 
approaches  have  been  empirical  and  have  dealt  only  with  cloud  fraction  (and  a  single  cloud  type  such  as 
cumulus;  e.g.,  Shenk  and  Salomonson  (1972),  Wielicki  and  Welch  (1986))  or  land  cover  classes  (e.g., 
Woodcock  and  Strahler  (1987),  Townshend  and  Justice  (1988)).  In  the  case  of  cloud  fraction,  real  and 
synthetic  data  containing  cloud  fields  were  degraded  in  resolution,  and  the  fractional  coverage  was 
observed  as  a  function  of  scale.  In  the  case  of  land  cover  classes,  the  variance  of  the  image  was  plotted 
as  a  function  of  measurement  scale  for  the  purpose  of  determining  the  optimal  resolution  for  monitoring. 
While  all  of  these  studies  are  useful,  no  concise  statement  of  the  relationship  between  fractional  coverage 
and  sensor  resolution  was  given,  so  that  the  results  are  difficult  to  generalize  to  other  parameters.  An 
analytical  solution  to  this  problem  is  currently  being  investigated,  but  here  we  present  results  from  an 
empirical  study. 

While  there  are  advantages  to  comparing  lead  statistics  derived  from  different  types  of  imagery, 
such  a  study  would  be  difficult  given  different  acquisition  times,  wavelength  ranges  of  the  various 
sensors,  and  geolocation  problems.  To  alleviate  these  problems,  comparisons  are  made  across  image 
scales  based  on  the  same  initial  image,  where  the  change  in  resolution  is  obtained  by  modeling  the  transfer 
function  between  the  initial  data  and  the  desired  resolution  and  then  subsampling.  A  spatial  filter  that 
estimates  the  point  spread  function  of  the  Landsat  sensor  is  applied  following  the  methodology  presented 
in  Justice  et  al.  (1989).  At  each  degradation  cycle,  Gaussian  random  noise  is  added  back  into  the  image 
to  reduce  the  smoothing  effects  of  the  filtering  operation.  Images  with  FOVs  of  160,  320,  640,  and  1280 
m  were  created  in  this  manner.  Data  are  Landsat  Multispectral  Scanner  (MSS)  band  4  (0.5-0.6  pm)  scenes 
of  the  Beaufort  Sea,  March  1988,  with  an  initial  field-of-view  (FOV)  of  <?0  m. 

Leads  statistics  are  retrieved  from  a  lead/not-Iead  binary  image.  To  obtain  the  binary  image  the 
fourth  order  trend  surface  is  removed  from  original  grey  scale  image  ( Eppler  and  Full ,  1992)  and  a 
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thresholding  procedure  is  applied  at  each  degradation  cycle.  Valid  lead  fragments  are  identified,  where 
"valid"  refers  to  a  linear  feature  for  which  a  meaningful  width  and  orientation  can  be  determined. 
Linearity  is  determined  through  correlation/regression  analysis.  Lead  widths  are  measured  perpendicular 
to  the  regression  line,  at  1  km  intervals,  and  the  slope  of  the  regression  line  is  the  measure  of  the  lead 
orientation. 

One  of  the  Landsat  scenes  analyzed  is  shown  in  Figure  3.  Image  size  is  80  x  80  km,  a  subset  of 
the  full  Landsat  scene.  The  corresponding  distribution  of  lead  widths  is  shown  in  Figure  4.  The 
disappearance  of  small  leads  due  to  a  lack  of  contrast  and  the  increase  in  the  relative  frequency  of  large 
leads  as  pixel  size  increases  can  readily  be  seen.  In  this  particular  image,  a  width  threshold  of 
approximately  250  m  divides  those  leads  which  would  disappear  in  the  320  and  640  m  degradations. 
However,  the  criteria  for  how  a  given  lead  will  "grow"  or  disappear  during  image  degradation  must 
include  a  measure  of  its  thermal  contrast  to  the  surrounding  ice.  For  example,  a  narrow,  open  lead  might 
grow  during  the  first  degradation,  accompanied  by  a  drop  in  contrast  due  to  smoothing,  and  then  disappear 
in  the  subsequent  scene.  A  narrow  refrozen  lead,  in  comparison,  might  disappear  during  the  first 
degradation.  This  change  in  lead  widths  as  a  function  of  pixel  size  affects  the  total  lead  area  as  shown 
in  Table  1  for  the  images  in  Figure  3,  where  lead  area  fraction  decreases  with  increasing  pixel  size. 
Orientations  of  leads  can  also  be  expected  to  change,  if  anisotropy  (i.e.,  a  preferred  orientation)  exists. 
An  illustration  of  this  is  shown  in  Figure  5  for  the  Landsat  image  in  Figure  3.  Results  from  other  Landsat 
scenes  show  similar  patterns. 

Given  that  very  small  features  will  generally  not  be  resolved,  the  issue  then  becomes  one 
concerning  the  possibility  of  using  the  partial  distribution  of  a  particular  statistic  (e.g.,  lead  widths) 
measured  at  low  resolution  to  estimate  the  complete  or  "true"  distribution  of  that  parameter.  For  example, 
assume  that  lead  widths,  x,  follow  a  negative  exponential  distribution  with  an  unknown  mean,  X.  From 
a  sampling  point  of  view  it  is  useful  to  treat  the  distribution  of  widths  as  discrete  and  address  the  number 
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Hi  of  leads  in  bin  i  that  have  widths  between  x,  and  xt+w  where  w  is  the  width  of  the  bin: 


Mv  -i,  A 

~Te 


(1) 


where  N  is  the  unknown  total  number  of  leads  in  the  spatial  area.  The  idea  is  that  n,  is  measured  for  a 
few  bins,  and  that  X  and  N  are  estimated.  To  accomplish  this,  (1)  is  rewritten  in  linear  form  as 


ln(n  )  =  \n(~)  ~  -i*,  (2) 

K  A 

Letting  a  =  ln(Nw/X)  and  b  =  X'1  and  solving  for  a  and  b  by  the  method  of  least  squares  with  the  observed 
data,  the  mean  of  the  distribution  and  the  total  number  of  leads  can  then  be  estimated. 

Experiments  with  this  model  show  it  to  be  very  sensitive  to  the  bin  width  and  the  number  of  bins 
in  which  leads  actually  occurred  in  the  low  resolution  imagery.  This  is  not  unexpected  considering  that 
the  entire  range  of  x  is  being  estimated  in  the  least  squares  model  by  observations  in  only  one  part  of  its 
entire  range  (such  extrapolation  is  not  recommended).  The  problem  can  be  alleviated  somewhat  by 
including  a  few  observations  of  small  lead  widths;  for  example,  measurements  from  a  few  SAR  images 
within  the  same  area. 

4.  SUMMARY  AND  CONCLUSIONS 

The  use  of  the  thermal  channels  of  the  AVHRR  for  the  detection  and  mapping  of  leads  has  been 
evaluated  for  varying  atmospheric  and  surface  conditions.  Our  approach  included  both  modeling  and 
empirical  studies,  discussed  with  respect  to  lead  width.  Concerning  atmospheric  effects,  the  detectability 
of  leads  is  dependent  on  the  degree  of  thermal  contrast  between  surface  types,  which  varies  with  surface 
and  atmospheric  conditions  as  well  as  sensor  response.  The  transfer  of  radiation  through  the  atmosphere 
is  modeled  using  data  from  Arctic  ice  islands  and  climatological  values  of  atmospheric  chemical 
constituents.  The  satellite  radiances  in  the  three  AVHRR  thermal  channels  are  simulated  for  a  variety  of 
surface  and  atmospheric  conditions:  a  range  of  satellite  view  angles,  four  ice  thicknesses  (three  within 
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leads),  varying  aerosol  and  ice  crystal  optical  depths,  and  thin  cirrus  clouds.  The  change  in  the  contrast 
between  the  lead  and  the  surrounding  ice  is  examined  as  a  function  of  optical  depth.  For  example,  at  ice 
crystal  optical  depths  of  more  than  0.3  under  mean  January  conditions  in  the  central  Arctic,  the  split- 
window  brightness  temperature  differences  between  2  m  and  5  cm  thick  ice  are  similar  to  the  ice 
temperature  variability.  Therefore,  thresholding  operations  will  not  be  able  to  distinguish  between  ice  and 
leads.  Due  to  the  wide  variety  of  atmospheric  and  surface  conditions  that  exist  in  the  Arctic,  additional 
experiments  in  this  area  must  be  performed  before  general  statements  about  the  determination  of  optimal 
thresholds  can  be  made. 

Concerning  measurement  scale,  Landsat  MSS  data  is  successively  degraded  from  80  m  to  1.2  km 
field-of-views,  using  the  modulation  transfer  function  for  the  sensor.  It  can  be  seen  that  small  leads 
disappear  in  the  coarser  resolution  data  and  large  leads  "grow".  The  total  lead  fraction  decreases  with 
resolution,  roughly  following  a  fractal  scaling  relationship.  Lead  orientation  distributions  also  change  with 
increasing  pixel  size  indicating  that,  in  the  imagery  examined,  small  leads  exhibit  different  orientations 
than  larger  leads.  The  possibility  of  using  lead  width  distributions  from  AVHRR  data  to  estimate  the 
"true"  lead  width  distribution  is  also  assessed,  where  the  large  leads  measured  by  satellite  are  first  adjusted 
for  their  artificial  growth,  and  then  their  distribution  is  used  to  predict  the  first  moments  of  the  true 
(hypothesized)  negative  exponential  distribution.  The  results  are  shown  to  be  very  sensitive  to  the  number 
and  range  of  lead  widths  measured,  but  the  accuracy  of  the  procedure  can  be  improved  dramatically  with 
a  few  samples  of  the  same  lead  network  at  a  higher  resolution;  e.g.,  from  SAR  data. 
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2407  (University  of  Washington  subcontract  #721566). 
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Table  1 .  Changes  in  lead  area  fraction  with  pixel  size  in  one  Landsat  scene. 


Pixel  Size  (m) 

- - - ■ - 

Lead  Area  Fraction 

Lead  Area  (kmJ) 

80 

0.030 

270.0 

160 

0.027 

235.8 

320 

0.023 

196.7 

640 

0.016 

124.9 
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Figure  1.  The  effect  of  a  layer  of  (a)  ice  crystal  precipitation  and  (b)  boundary  layer  haze  on  the  three 
AVHRR  thermal  channels  (3,  4,  5)  at  a  satellite  view  angle  of  50°.  Optical  depths  (0.55  pm)  are  varied 
from  0  to  10.  The  tops  of  the  ice  crystal  and  haze  layers  are  near  the  top  of  the  inversion  at  a  temperature 
of  approximately  248  K.  The  brightness  temperature  change  in  each  channel  over  each  of  the  four  surface 
types  (ice  thicknesses:  0,  5,  15,  200  cm  from  top  to  bottom)  can  readily  be  seen. 

Figure  2.  Normalized  contrast  between  an  open  water  lead  and  ice  of  two  thicknesses  for  (a)  ice  crystal 
precipitation  and  (b)  boundary  layer  haze.  Conditions  are  mean  January  in  the  central  Arctic. 

Figure  3.  Landsat  MSS  band  4  scene  of  the  ice  pack  north  of  Alaska  in  March  1988.  Area  covered  in 
the  upper  left  image  is  approximately  (80  km)2.  The  remaining  are  degraded  images  with  pixel  sizes  of 
160,  320,  and  640  m. 

Figure  4.  Lead  width  distributions  for  the  images  in  Figure  3.  Widths  are  grouped  in  100  m  bins. 

Figure  5.  Lead  orientations  for  the  degraded  Landsat  series  shown  in  Figure  3.  Orientation  is  the  angle 
that  a  lead  makes  with  the  horizontal  axis,  measured  counterclockwise. 
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